AD-A285  315 


L/ 

AD 


LONG  TERM  MODELLING 
OF 

PERMAFROST  DYNAMICS 


4^-n-  cTK-Di 

T>TiC 


Final  Technical  Report 


by 

K. Fawcett  and  M.G. Anderson 


July  1994 


European  Research  Office 
U.S.  Corps  of  Engineers 
London  England 


V  . 


CONTRACT  NUMBER  DAJA  45-92-C-001 1 

Professor  M.G. Anderson 


Approved  for  Public  Release  ;  Distribution  Unlimited 


Best 

Available 

Copy 


The  PERIMOVE  model  is  capable  of  simulating  moisture  and  heat  transfers  within  soil  profiles. 
In  addition,  the  simulation  of  these  processes  formed  the  basis  for  the  prediction  of  the  volume 
changes  that  moisture  redistribution  and  freezing  create  within  soil  profiles  and  the  resulting 
movement  of  material  and  loss  of  soil  strength  during  freezing  and  thawing  cycles.  In  order  to 
aid  the  development  and  validation  of  the  model  a  suite  of  data  sources  have  been  identified  and 
assembled.  These  have  enabled  model  testing  at  each  stage  of  development.  Areas  of  further 
work  are  suggested. 

1  ORIGINAL  OBJECTIVES 

ORIGINAL  OBJECTIVE; 

To  develop  a  2-dimensional,  physically  based  model  for  mass  and  heat  transfer,  driving  a 
creep  process  model  for  cold  regions  slopes. 

Harris  et  al  (1993)  describe  some  results  of  a  NERC  funded  project  involving  a  large  scale 
laboratory  experiment  into  frost  creep  and  gelifluction  processes.  These  processes  can  be 
identified  as  distinct  where  creep  describes  the  downslope  movement  of  material  due  to 
gravitational  resettlement  of  soil  during  thawing  as  a  result  of  heave  occurring  during  a  freeze 
cycle  and  gelifluction  describes  the  downslope  deformation  of  a  saturated  mass  of  soil  overlying 
a  frozen  soil  layer  during  the  thaw  cycle.  Figure  1  shows  their  results  for  four  different  soil  types 
and  the  relative  importance  of  creep  and  gelifluction  processes  in  controlling  downslope 
movement  for  the  different  soil  types  is  apparent.  Because  of  the  significance  of  gelifluction  in 
producing  downslope  movement  of  material  it  was  decided  that  the  original  objective  should  be 
altered  to  include  gelifluction  processes  in  addition  to  creep  processes. 

The  term  solifluction  has  replaced  the  term  creep  in  the  original  objective  as  the  former  term 
embraces  both  frost  creep  and  gelifluction  movements  as  defined  in  Harris  (1987).  Collaboration 
of  this  project  with  that  of  Harris  has  taken  place  and  is  described  in  more  detail  in  later  sections. 

An  additional  change  in  the  original  objective  is  that  the  model  should  strictly  be  termed  a 
pseudo  2-dimensional  model.  Although  the  basic  heat  and  moisture  transfer  model  could  operate 
in  2-dimensions  the  increase  in  complexity  resulting  from  the  incorporation  of  volume  change 
make  a  fully  2-dimensional  scheme  beyond  the  means  of  this  project.  Ideally  the  model  would 
operate  within  a  continuously  deforming  2-dimensional  mesh  structure.  Possibilities  for  the 
extension  to  this  stracture  are  discussed  in  section  6.  PERIMOVE  in  its  current  form  however 
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is  capaole  of  simulating  the  2-dimerisional  transfers  of  material  which  occur  over  freeze/thaw 
cycles. 


2  DEVELOPMENT  OF  METHODS  AND  TECHNIQUES 

Three  areas  of  development  are  outlined  below.  The  first  section  outlines  the  structure  and 
equation  base  of  the  PERIMOVE  model.  The  second  section  outlines  the  data  preprocessing 
scheme  and  the  third  section  describes  the  data  sources  that  have  been  collated  to  supply  data 
which  have  contributed  to  the  development  and  the  validation  of  the  model.  Figure  2 
summarises  the  development  and  validation  route  for  the  model  and  also  suggests  to  further  areas 
of  development  for  the  model. 

i)  Description  of  processes,  equations  and  spatial  discretisation  approach 
selected  for  the  PERIMOVE  model 

Spatial  Discretisation 

The  slope  profile  is  represented  as  a  series  of  columns  of  equal  width  and  breadth.  Each  column 
is  further  subdivided  to  create  cells  which  are  then  treated  as  homogeneous  units  within  the 
model.  The  computation  point  is  at  the  centre  of  each  cell  and  all  thermal  and  hydraulic 
properties  are  assigned  to  that  point.  These  block-centred  computation  points  represent  the 
mathematical  nodes  for  the  finite  difference  mass-energy  calculations  performed  by  the  model. 
Mass  and  energy  transfers  take  place  on  the  basis  of  the  thermal  and  hydrological  gradients 
which  are  calculated  between  cells  and  phase  change  processes  are  governed  by  the  within  cell 
conditions. 

A  maximum  of  two  hydrologically  distinct  layers  may  be  identified  (as  defined  by  suction- 
moisture  curves  and  characteristic  freezing  curves).  The  boundary  between  them  is  determined 
by  specification  of  the  lowest  cell  in  the  upper  layer  so  that  the  rest  of  the  column  comprises  the 
second  layer. 

Each  slope  section  is  modelled  in  isolation  therefore  three-dimensional  attributes  such  as  aspect, 
convergent  and  divergent  slope  forms  are  not  included  in  the  present  model. 

Process  and  Equation  Base 

Three  major  sets  of  equations  can  be  identified  which  describe  the  three  major  groups  of 
processes  which  need  to  be  included  in  a  modelling  scheme.  These  are: 
processes  of  heat  flow  in  soils 

processes  of  phase  change  in  freezing  and  thawing  soils 
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MODEL  DEVELOPMENT 

-  spatial  discretization 

-  processes  of  heat  exchange 

-  processes  of  phase  change 

-  volume  change 

simulate  with  fully 
flexible  mesh  structure 


SECTION  2(i) 


INPUT  DATA  REQUIREMENTS  ! 

-  values  from  literature 

-  CRREL  data  sources  j 

-  data  preprocessor  (STEPS)  ! 


SECTION 

2(u) 


Collate  data  sources 
for  validation 
SECTION  2(ii) 


determine  validation 
data  required 


MODEL  VALIDATION  ! 

-  model  behaviour:  SENSITIVITY  | 

-  heat  transfer  and  soil 

heave:  FERE 

-  heat,  moisture  transfer 

and  soil  displacement: 


SECTION 


CAEN 


suggest  developments 


Fully  integrated  2-dimensional  solifluction 
model  operating  within  a  deforming  finite 
element  mesh 


CAEN  -  fully  instrumented 

slope  data  —  1994/95 
Field  validation  data 


Figure  2;  Summary  diagram  of  stages  in  model  development  and  application.  Solid  line 
boxes  indicate  the  scope  of  this  project  and  dotted  lines  indicate  areas  for  further  research. 
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processes  of  water  transfer  through  the  profile. 

Although  these  processes  all  interact  it  is  possible  to  consider  them  independently  in  order  to 
outline  the  major  equations  describing  each  of  the  processes. 

Processes  of  heat  flow  within  soil 

Heat  flow  within  soils  takes  place  by  a  variety  of  processes.  In  periglacial  situations  the  two 
major  processes  responsible  for  heat  transfer  are  conduction  through  the  soil/water/ice  medium 
and  transfer  of  latent  heat  by  liquid  flow.  Additional  processes  such  as  heat  flow  by  radiation, 
heat  generated  by  mechanical  work  of  the  particles  and  sensible  heat  convection  are  generally 
considered  to  be  of  minor  importance  and  are  currently  excluded  from  the  PERIMOVE  model. 

Thermal  Conductivity  is  a  measure  of  the  quantity  of  heat  that  will  flow  through  a  unit  area  of  a 
substance  of  unit  thickness  in  unit  time  under  a  temperature  gradient.  It  is  usually  measured  in 
Wm-1  K-l. 

In  situations  where  ground  temperatures  fall  below  0°C  the  thermal  conductivity  of  soils  cannot 
be  assumed  constant  as  freezing  soils  have  a  temperature  dependent  ice  content  and  the  thermal 
conductivity  of  ice  is  four  tmes  that  of  water.  In  soils  undergoing  freeze/thaw  the  ice/water 
ratios  change  rapidly.  In  periglacial  situations  where  soil  temperatures  are  frequently  in  the 
range  0  to  -3°C,  failure  to  account  for  the  temperature  dependence  of  thermal  conductivity  may 
lead  to  unacceptable  errors  in  the  modelling  stage  so  temperature  dependant  relationships  are 
incorporated  in  the  PERINIOVE  model. 

DeVries  (1952)  developed  the  equation  for  calculation  of  thermal  conductivity  in  a  multi- 
component  system: 

Kt  =  .  I  Xi  kt:  F: 

n 

.  E  Xi  Fi 
1=0  ‘  ^ 

(eq.  1) 

Kt  =  thermal  conductivity  of  the  system  (W  m'^  K'^) 

n  =  number  of  different  kinds  of  particles  in  the  continuous  media 

Xj  =  volume  fraction  of  ith  kind  of  particle  (%) 

k[j  =  thermal  conductivity  of  ith  kind  of  particles  (W  m'  ^  K" 

Fj  =  ratio  of  average  temperature  in  the  ith  kind  of  particles  to  average 
temperature  in  the  continuous  media 

this  means  that  for  a  water,  ice  and  soil  system  the  equation  is  (by  Penner,  1970): 
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Kt  +  XjktjFj  +  XjktgFjj 


(eq.  2) 


subscripts; 


w  =  water,  (F^=l) 

i  =  ice 

s  =  soil 


PERIMOVE  uses  the  above  formula  in  which  both  ice  and  water  thermal  conductivities  are 
considered  temperature  dependent  below  0°C.  The  relationships  of  thermal  conductivity  with 
temperature  are  input  to  the  model  as  a  series  of  data  pairs  of  conductivity  and  temperature  and 
interpolation  between  data  points  is  made  to  'look-up'  conductivity  values  for  the  current 
temperature  of  the  particular  computation  cell. 

The  above  equation  does  not  include  latent  heat  transfer  effects  and  as  stated  above  significant 
errors  can  occur  in  predictions  if  these  are  omitted.  Latent  heat  is  released  during  freezing  and 
absorbed  during  thawing  and  means  that  a  partially  frozen  soil  behaves  as  if  it  has  a  greatly 
increased  heat  capacity.  Latent  heat  is  considered  in  the  following  section. 

Once  a  therm.al  conductivity  has  been  established  it  can  be  used  to  describe  transfer  of  heat  by 
conduction  by  application  of  Fourier's  Law; 

h^  =  -Kt^  5TA 


where; 


h^  =  rate  of  energy  flow  in  x  direction  (W) 

Kt  =  thermal  conductivity  in  specified  direction  (W  m'^  K'^) 
T  =  temperature  (K) 

A  =  area  (m^) 


Fourier's  law  states  that  the  rate  of  energy  flow  in  the  direction  of  temperature  drop  is 
proportional  to  the  area  through  which  the  energy  flows  and  the  temperature  gradient. 


Laplace's  equation  combines  the  heat  conduction  theory  with  the  continuity  equation  to  generate 
a  second  order  partial  differential  equation; 


5  /'-Kt5h"j  +  5  r-Kt5h]  +  5  f-Kt  6h '|  =  0 

5x  5x J  5y  V  5y  j  5z  5z  ^ 

(eq.  3) 

where;  Kt  =  Kt  (x.y.z) 

If  the  region  considered  is  homogeneous  and  isotropic  then  Kt  is  independent  of  x,y  and  z,  so  that 
the  above  equation  simplifies  to: 


•Kt5hl  H- 


-Kt  6h 
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5“  h  +  §2  h  +  5^  h  =  0 
5x2  5y2  5^2 

(eq.  4) 


Laplace's  equation  therefore  describes  heat  conduction  under  steady  state  conditions.  Fourier's 
law  in  combination  with  the  continuity  equation  is  solved  in  one  dimension  within  PERIMOVE 
so  that  temperature  change  through  time  is  described  by: 


Cv  5T 
5t 


-5  ["-Kt  5T' 
5x  ^  5x  j 


where:  Cy  =  volumetric  heat  capacity  of  the  soil  medium 


(eq.5) 


and:  volumetric  Heat  Capacity  is  the  amount  of  heat  required  to  change  the  temperature  of  a 
unit  volume  of  a  substance  by  one  degree.  The  heat  capacity  of  a  mixture  is  equal  to  the  sum  of 
the  heat  capacities  of  its  components. 


This  heat  exchange  with  the  profile  is  achieved  by  specification  of  Neuman  boundary  conditions 
in  which  the  transfers  of  energy  across  the  surfaces  bounding  the  problem  are  known. 


The  thermal  response  of  a  soil  requires  the  calculation  of  volumetric  heat  capacity.  In  unfrozen 
soils  the  volumetric  heat  capacity  is  dependent  on  dry  bulk  density  of  soil,  moisture  content  and 
specific  heat  capacity  of  the  constituent  mineral  matter,  (Jumikis,  1966).  The  relationship  can  be 
expressed  as: 


Cy  =  yd  (  Cs  +  Cw  W/100) 


where: 

Cy  =  volumetric  heat  capacity  (J  m"^ 
yd  =  dry  bulk  density  of  soil  (g  m'^) 

Cs  =  specific  heat  capacity  of  dry  mineral  matter  (J  g'^  K'^) 
Cw=  specific  heat  capacity  of  water  (4.18J  g'^  K'^) 

W  =  water  content  (%) 


As  soil  freezes  the  volumetric  heat  capacity  changes  as  relative  volumes  of  ice  and  water  change. 
The  equation  for  partially  frozen  soil  is  expressed  by: 

Cyf  =  ydf  (  Cs  +  (Cw  W/100)  +  (Ci  I/lOO) ) 

where; 

ydf  =  frozen  bulk  density 

Ci  =  specific  heat  capacity  of  ice  (J  g'^  K'^)  pure  ice  =  2.1  J  g‘^  K"^ 

I  =  ice  content  (%) 
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This  assumes  that  the  air  phase  can  be  ignored  and  the  specific  heats  are  assumed  to  be 
independent  of  temperature.  Bi-modal  values  of  soil  bulk  density  where  one  value  is  used  for 
frozen  conditions  and  another  for  unfrozen  conditions  are  employed  m  PERIMOVE.  This 
represents  a  simplification  of  reality  as  soil  structure  changes  during  repeated  freezing  and 
thawing  events  producing  corresponding  changes  in  bulk  density  values.  It  may  well  be 
necessary  to  incorporate  this  continual  change  in  bulk  density  at  a  later  stage  especially  if  the 
model  is  particularly  sensitive  to  this. 


Processes  of  Phase  Change  in  Freezing  and  Thawing  Soils 

Phase  change  of  water  to  ice  and  ice  to  water  takes  place  over  a  range  of  temperatures  in  soils. 
As  an  example,  fine-grained  silt/clay  soils  can  still  contain  a  liquid  moisture  content  at 
temperatures  as  low  as  -20°C.  Therefore  it  is  not  sufficient  to  consider  the  freezing  process 
occurring  at  a  single  interface  (i.e.  the  0°C  isotherm).  Phase  changes  in  cooling  and  warming 
soils  are  determined  by  heat  flux  and  water  content  of  the  cell. 

Phase  change  in  freezing  soils 

The  heat  budget  for  a  uniform  volume  of  soil  can  be  described  by: 

AQ  =  La5e,  AT  -  C^5AT 

5T  (eq.  6) 

where; 

La  59/5T.AT  is  the  energy  released  due  to  phase  change  of  moisture  to  ice, 
termed  the  hydrological  component. 

CySAT  is  the  energy  release  due  to  thermal  capacity  of  soil  medium, 

termed  the  thermal  component. 

AQ=  change  in  energy  of  the  system 
La  =  latent  heat  of  fusion  (J  g'^) 

0  =  volumetric  soil  moisture  content  (%) 

Cy  =  volumetric  heat  capacity  (J  m'^  K'^) 

V  =  volume  (of  cell)  (m'^) 

T  =  temperature  (K) 

Use  is  made  of  the  characteristic  freezing  curve  for  a  soil  which  is  an  experimentally  determined 
curve  which  describes  the  percentage  of  liquid  water  which  remains  unfrozen  at  a  given 
temperature  below  0°C 

The  characteristic  freezing  curve  for  the  specified  soil,  the  soil  liquid  water  content  and  the  soil 
temperature  are  used  to  determine  the  quantity  of  ice  formed  for  a  given  drop  in  temperature. 
The  operation  of  the  curve  in  the  PERIMOVE  model  is  described  in  figure  3. 


Characteristic  Freezing  Cup;e 


I  The  characteristic  freezing  curve  for  a  soil  defines  the  relationship  between  liquid 
i  water  content  and  temperatures  below  273. 15K  (O’C).  The  cur/e  can  be 
'  determined  experimentally  and  an  example  of  a  characteristic  freezing  curve  is 
;  given  below. 

Unfrozen  moisture  content  (”0) 


6 1 , 


:  if  the  soil  is  initially  saturated  (Gg)  then  as  the  soil  temperature  begins  to  drop 
;  below  O’C  the  liquid  water  content  will  fall  such  that  at  temperature  T-|  the  liquid 
;  water  content  will  be  Gi  and  at  T3  the  corresponding  liquid  water  content  will  be 
I  83.  The  fall  in  liquid  water  content  with  temperature  is  equivalent  to  the  water 
;  equivalent  volume  of  ice  formed. 

!  If,  however,  the  soil  is  unsaturated  at  the  start  of  the  freezing  prcccsc,  for  example 

( 

;  02  then  ice  formation  will  not  be  initiated  until  the  soil  temperature  falls  below  T2 
(the  temperature  corresponding  to  a  liquid  water  content  of  82).  Untii  this  point  is 
reached  the  liquid  water  content  will  be  super-cooled.  Once  temperature  T2  is 
reached  the  falling  liquid  water  content  will  follow  the  pattern  of  the  characteristic 
curve. 

This  curve  is  used  in  subroutine  FREEZE  to  determine  a  liquid  water  content 
given  the  soil  temperature  (as  calculated  in  HEATRANS). 


Figure  3:  Operation  of  the  characteristic  freezing  curve  within  the  PERIMOVE 

model. 
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On  freezing  thr. .  possible  thermo-hydrodynamic  conditions  are  considered. 

1)  La  1'*^  >  Cy  AT  AQ  is  positive  and  isothermal  phase  change  occurs 

5T 

The  energy  generated  from  the  phase  change  of  water  into  ice  is  greater  than  the  temperature 
decline  controlled  by  ambient  temperature  change,  heat  conduction  through  the  soil  and  heat 
capacity  of  the  soil  (thermal  component).  The  soil  temperature  remains  unchanged  as  the 
temperature  decline  due  to  soil  heat  capacity  and  temperature  gradient  is  compensated  for  by 
energy  release  from  phase  change.  Excess  AQ  is  assumed  to  be  lost  from  the  system. 

2)  La  50  AT  =  Cy  AT  AQ  is  zero 

5T 

The  energy  losses  due  to  the  thermal  component  are  equivalent  to  the  energy  gain  due  to  water 
phase  change.  The  soil  temperature  remains  constant  as  the  gains/losses  are  compensatory. 

3)  La  50  AT  <  Cy  AT  AQ  is  negative,  non-isothermal  phase  change  occurs 

5T 

The  energy  loss  from  the  thermal  component  is  greater  than  can  be  compensated  for  by  release  of 
energy  which  occurs  as  the  soil  cools  below  freezing.  There  is  a  drop  in  soil  temperature  as 
phase  change  cannot  release  sufficient  energy.  AQ  is  negative  and  no  energy  is  lost  from  the  soil 
system  as  all  hydrological  latent  heat  is  used  in  compensating  for  thermal  losses. 


Thawing  Soils; 

The  logic  of  phase  change  for  thawing  soils  is  based  on  the  assumption  that  energy  input  to  a  cell 
is  partitioned  between: 

(1)  temperature  change  -  the  thermal  component 

(2)  phase  change  of  ice  to  water  -  the  hydrological  component 

A  thaw  temperature  is  specified  (normally  273. 15K)  and  all  energy  input  into  a  cell  is  used  to 
raise  the  cell  temperature  until  the  thaw  temperature  is  reached.  Once  the  thaw  temperature  is 
reached  the  energy  transferred  into  a  cell  is  used  to  satisfy  the  energy  requirements  of  the  phase 
change  from  ice  to  liquid  water.  Once  the  cell  is  ice  free  the  incoming  energy  is  again  used  to 
raise  the  cell  temperature.  These  can  be  expressed  as  three  conditions. 

1)  if  the  cell  temperature  is  less  than  the  set  thaw  temperature  then  all  the  energy  entering  the 
cell  is  used  to  raise  the  cell  temperature. 


AT=  AQ/  CySAT 


2)  if  cell  temperature  is  at  thaw  temperature  and  ice  is  present  then  the  energy  entering  the  cell 
IS  used  to  satisfy  the  latent  heat  requirement  of  melting  ice.  The  incoming  energy  determines  the 
volume  of  ice  which  can  be  melted.  The  temperature  remains  at  the  thaw  temperature  until  all 
the  cell  ice  is  melted. 


A0=  AQ  /  La59  .  AT 

6T 


3  i  if  the  cell  temperature  is  at  or  greater  than  thaw  temperature  and  there  is  no  ice  present  all 
the  incoming  energy  is  used  to  raise  the  cell  temperature. 

AT  ^  AQ  /  CySAT 


Processes  of  Water  Transfer  Through  the  Profile 

In  the  current  version  of  PERCVIOVE  only  vertical  movement  of  water  is  incorporated. 

Darcy's  law  provides  the  basic  equation  for  water  flow  through  soils  and  for  vertical  movement 
the  equation  can  be  written  as: 


qy=  -K  6H  =-K5  (Hp-y) 

5y  5y 

where; 

qy  =  water  flux 
H  =  total  hydraulic  head 
Hp  =  pressure  head  (m) 

y  =  vertical  distance  component  (in  PERIMOVE  scheme) 
K  =  hydraulic  conductivity  (m  s' ') 


(eq.  7) 


In  unsaturated  soils  Hp  is  negative  and  can  be  expressed  as  a  suction  head; 

qy  =  K  5v/  +  K 

5y 


Combining  the  above  equation  with  the  continuity  equation  yields: 


(eq.  8) 


Ai  soil  moisture  content  and  suction  head  are  uniquely  related  the  left  hand  side  of  the  equation 
can  be  written  as: 

5t  5v)/  5t 


Substituting  this  relationship 

C5m/  =  5  f'  K6m/' 

5t  5y[  5y  ^ 


into  the  above  equation  produces: 

+  5K 

.5y 


(eq.  9) 


where:  C  is  the  specific  water  capacity  describing  the  change  in  water  content  of  a  unit  volume 
of  soil  per  unit  change  of  matric  suction. 

Equations  8  and  9  are  forms  of  the  Richards  equation.  Conventionally  these  are  developed  for 
flow  in  isothermal  saturated  and  unsaturated  flow  within  ice  free  soil  water  systems. 

Extension  of  this  equation  to  situations  where  ice  formation  occurs  within  the  soil  is  possible 
because  of  the  similarity  between  air/water  systems  and  ice/water  systems.  As  temperatures  fall 
to  the  point  of  ice  nucleation  a  single  crystal  forms  in  the  centre  of  a  water  filled  pore.  If  this 
crystal  is  regarded  as  analogous  to  an  air  bubble  then  the  flow  solution  can  be  assumed  similar  to 
that  for  unsaturated  flow.  The  growth  of  an  ice  crystal  can  be  compared  to  the  growth  of  an  air 
bubble  in  ursaturated  soils  and  relationships  developed  for  drying  soils  can  be  used  to  describe 
the  hydraulic  properties  of  freezing  soils  (Burt  and  Williams,  1976). 

As  ice  begins  to  form  in  the  upper  soil  layers  the  liquid  moisture  content  of  these  layers  falls. 
This  results  in  significant  soil  water  suctions  developing  in  these  cells  and  allows  for  water  to  be 
drawn  up  into  the  freezing  zone  and  fuel  the  development  of  ice  quantities  which  can  cause 
profile  expansion.  Soil  conductivity,  moisture  availability  and  rate  of  freezing  control  the  extent 
to  which  water  will  be  drawn  up  the  profile. 

The  hydraulic  conductivity  curve  with  moisture  content  is  derived  using  the  Millington-Quirk 
procedure  within  the  PERIMOVE  model. 


Volume  change  and  mass  movement  in  the  PERIMOVE  model 

The  incorporation  of  volume  expansion  and  the  movement  of  water  through  the  profile  is 
essential  if  volume  expansion  during  freezing  and  subsequent  supersaturation  of  the  soil  profile 
and  loss  of  soil  strength  following  thawing  are  to  be  modelled.  It  is  known  that  volume 
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expansion  of  soil  undergoing  freezing  is  not  due  solely  to  the  volume  change  that  occurs  when 
liquid  water  becomes  ice  within  the  soil  pore  system.  Continued  migration  of  water  towards  the 
freezing  front  as  liquid  water  is  converted  to  ice  results  in  the  formation  of  segregation  ice  lenses 
and  additional  expansion.  On  thawing  this  excess  water  is  often  prevented  from  draining  away 
due  to  the  frozen  state  of  lower  profile  layers.  The  supersaturation  of  the  soil  material  that  results 
creates  conditions  which  can  lead  to  solifluction  because  of  the  reduction  in  soil  strength. 

For  inclusion  in  the  PERIMOVE  model  it  is  assumed  that  the  change  from  liquid  water  to  ice  is 
accompanied  by  a  9%  increase  in  volume.  Within  each  cell  the  ice  formation  proceeds  until  the 
sum  of  the  liquid  water  and  ice  contents  exceeds  the  porosity  of  the  soil.  When  this  situation  is 
reached  the  new  volume  of  the  cell  and  the  change  in  cell  depth  (all  expansion  is  assumed  to  take 
place  in  the  vertical  dimension)  is  calculated.  Following  this  the  new  volumetric  percentages  of 
soil  particles,  liquid  water  and  ice  are  calculated. 

The  new  values  are  then  transferred  to  the  mass  transfer  calculation  section  of  the  program  and 
transfers  are  made  based  on  the  hydraulic  gradients  between  cells. 

Previously  the  movement  of  water  between  cells  was  limited  by  the  cell  porosity  as  water  could 
not  move  into  a  cell  which  had  been  'filled'.  This  check  has  been  removed  to  allow  continued 
upward  migration  of  water  to  fuel  freezing  and  volume  expansion. 

The  volume  expansion  of  a  cell  also  effects  the  operation  of  the  suction-moisture  curve  and  the 
characteristic  freezing  curves. 

Characteristic  freezing  curve 

In  normal  operation  the  characteristic  freezing  curve  is  used  to  determine  the  liquid  water  content 
of  a  cell  at  a  given  temperature  and  determine  if  water  is  available  in  the  cell  for  the  freezing 
process  to  proceed.  With  volume  expansion  the  assumption  is  made  that  the  liquid  water  at  a 
given  temperature  is  an  absolute  volume  calculated  on  the  basis  of  the  original  cell  volume. 
When  volume  change  occurs  the  characteristic  freezing  curve  is  used  to  look  up  the  water 
content  for  the  given  temperature  and  a  value  for  volumetric  water  content  value  is  returned.  In 
the  program  the  value  is  then  transformed  to  a  value  in  m^  for  the  freezing  calculations.  If  the 
volume  change  is  not  accounted  for  in  this  calculation  then  the  volume  in  m^  will  be 
overestimated.  To  avoid  this  the  value  of  liquid  water  obtained  form  the  characteristic  freezing 
curve  is  divided  by  the  total  volume  change  that  has  occurred  over  the  simulation. 

Suction  moisture  curve 

The  suction  moisture  curve  is  used  to  'look-up'  values  of  m  suction  for  a  given  soil  moisture 
content.  The  curve  is  derived  for  liquid  moisture  contents  so  as  ice  forms  in  a  cell  liquid 
moisture  contents  fall  and  suctions  increase. 

As  volume  expansion  occurs  the  same  m3  volume  of  water  will  be  equivalent  to  a  reduced 
volumetric  water  content  and  return  a  higher  suction  value  from  the  suction-moisture  curve. 
However,  this  interpretation  may  be  satisfactory  as  an  equivalent  absolute  volume  in  a  larger  cell 
volume  will  lead  to  a  reduced  connectivity  in  pore  space  water  and  thus  a  higher  suction.  This  is 
the  approach  that  is  being  taken  in  the  current  version  of  the  PERIMOVE  model. 
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Conductivity 

The  Millington-Quirk  relationship  is  used  to  derive  a  conductivity-moisture  curve  from  the 
suction-moisture  curve.  The  conductivity-moisture  curve  is  then  used  to  'look-up'  the  value  of 
conductivity  for  a  given  volumetric  moisture  content  in  the  mass  transfer  section  of  the  program. 
If  volume  expansion  has  occurred  the  volumetric  liquid  water  content  will  be  lower  for  a  given 
absolute  volume  of  water  and  a  lower  conductivity  will  be  returned  than  for  the  same  liquid 
volume  pre-expansion.  This  reflects  the  decrease  in  connectivity  that  will  result  when  a  cell 
expands  so  the  values  extracted  from  the  curve  are  used  directly  in  the  model. 

Summary 

The  above  sections  have  described  the  model  structure  and  equations.  Figure  4  illustrates  the 
sequence  in  which  the  processes  described  above  are  modelled  within  each  time  step  of  the 
PERIMOVE  model. 

The  other  factors  to  be  considered  in  this  section  are  the  data  sources  that  have  been  assembled  to 
help  construct  and  validate  the  model. 

ii)  Soil  type  driven  data  preparation  scheme 

.A.n  important  part  of  the  model  development  plan  was  that  the  mo  !  would  successfully 
reproduce  the  differences  in  moisture  redistribution  and  freezing  of  different  soil  types.  To 
facilitate  the  derivation  of  data  files  for  this  the  Soil  Type  Essential  Property  Selector  (STEPS) 
program  was  written.  This  is  based  on  the  model  written  by  Cochrane  and  Anderson  (1988)  and 
generates  soil  property  files  when  given  a  soil  class  from  the  standard  US  DA  soil  triangle.  The 
original  model  and  database  have  been  expanded  to  include  selection  of  soil  thermal  properties 
and  characteristic  freezing  curves.  This  provides  a  quick  and  user  friendly  way  of  selecting  soil 
properties  for  the  PERIMOVE  model  and  introduces  the  option  of  wide  applicability  without 
elaborate  testing  of  soils. 

iii)  Validation  data  sources 

Two  major  data  sources  have  been  used  in  the  model  development  and  validation.  They  are 
derived  from  two  areas  of  collaboration  that  have  been  established  during  the  project. 

A)  Collaboration  with  the  Cold  Regions  Research  and  Engineering  Laboratory  in 
Hanover,  New  Hampshire.  A  range  of  experiments  conducted  there  have  been  useful  in 
development  and  testing  for  the  PERIMOVE  model. 

CRREL  has  conducted  many  laboratory  tests  to  derive  characteristic  freezing  curves.  In  a  visit 
made  there  the  various  sample  results  were  collated  and  soil  type  information  matched  with  the 
relevant  curves  to  augment  the  data  preprocessing  program  STEPS  described  above. 
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Figure  4;  Sequence  of  operation  of  the  subroutines  within  the  PERIMOVE  model. 
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Shoop  (1989)  describes  the  results  of  an  experiment  performed  within  the  Frost  Effects  Research 
Facility  (FERF)  at  CRREL.  The  FERF  building  is  an  environmentally  controlled  building  and 
contains  a  series  of  test  cells  measuring  around  20ft  square  and  12  ft  deep.  Soil  freezing  is 
instigated  by  placing  freezing  panels  on  the  soil  surface  and  thawing  is  produced  by  either 
heating  the  building  or  opening  the  building  to  the  ambient  air  temperature.  The  data  which  was 
made  available  was  collected  during  a  freeze/thaw  cycle  for  a  silty  sand.  Profile  temperatures 
were  monitored  and  measurement  of  total  soil  heave  made  at  the  end  of  the  freeze  cycle. 
Moisture  content  was  measured  using  tensiometers  (although  these  cease  working  once  freezing 
commences)  and  gravimetric  moisture  determinations  were  made  on  samples  of  rofile 
removed  at  the  end  of  the  freeze  cycle.  As  such  the  data  provides  a  good  source  for  ..ludating 
temperature  and  moisture  movernent  within  the  profile. 

An  additional  data  source  that  was  used  more  in  the  early  development  of  the  model  was  data 
collected  from  some  plot  experiments  conducted  at  CRREL.  Two  outdoor  plots  -  one  containing 
a  sandy  soil  type  and  the  other  a  loam  soil  type  were  monitored  over  two  winter  seasons.  The 
data  consists  predominantly  of  temperature  data  and  this  was  used  in  the  early  development  of 
the  model  before  other  data  sources  were  available.  The  second  season  has  some  moisture  data 
and  may  be  used  at  a  later  date. 

B)  The  second  source  of  data  has  come  from  collaboration  with  Harris  et  al  (1993)  who 
conducted  a  large  scale  experiment  using  a  constructed  slope  within  a  refrigerated  container 
(some  results  of  which  are  described  in  section  1)  using  different  soil  types  at  the  C.N.R.S. 
Centre  de  Geomorphologie  in  Caen,  France.  The  results  from  this  study  include  temperature 
profiles  and  heave  measurements  for  a  suite  of  freeze/thaw  cycles.  The  heave  measurements  can 
be  used  to  indicate  the  moisture  redistribution  occurring  during  freezing  within  different  soil 
types  and  the  temperature  profiles  to  compare  with  model  simulated  temperatures.  In  addition 
measurements  were  made  of  the  relative  amounts  of  downslope  movement  of  material  that  could 
be  attributed  to  creep  and  gelifluction  processes  for  each  of  the  soils  types.  This  gives  an 
extremely  good  data  set  against  which  soil  strength,  soil  supersaturation  on  thawing  and  soil 
movement  can  be  calibrated  within  the  model. 

3  PROGRESS  IN  RELATION  TO  SCHEDULE  OF  WORK 

At  the  current  stage  of  development  the  model  can  be  used  to  investigate  the  nature  and  extent  of 
2'dimensional  movement  of  material  downslope  although  as  stated  in  section  (1)  the  current 
model  is  pseudo  2-dimensional  in  structure  because  of  the  intense  computational  difficulties  in 
defining  a  2-dimensional  deforming  mesh  structure.  The  continued  development  of  finite 
element  techniques  will  mean  that  it  should  be  possible  in  the  future  to  have  a  fully  integrated  2- 
dimensional  simulation  capacity. 

The  collation  of  the  range  of  data  sources  for  model  validation  and  testing  has  been  an  important 
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part  of  the  project  and  a  larger  part  than  originally  anticipated  however,  the  collection  and  use  of 
these  data  sources  has  meant  the  model  development  and  validation  has  been  more  rigourous 
than  originally  planned. 


4  SUMMARY  OF  RESULTS  AND  DATA  TO  DATE 

The  results  summary  is  presented  in  three  sections  and  each  section  highlights  a  particular  area  of 
results  from  the  model.  The  results  follow  a  progression  through  the  validation  and  application 
development  of  the  model.  The  first  section  looks  at  the  general  performance  and  behaviour  of 
the  model  by  summarising  the  results  of  the  sensitivity  analysis  of  model  parameter  values.  The 
second  and  third  sections  concentrate  on  the  application  of  the  model  to  the  two  test  data  sources 
looking  first  at  the  application  to  the  FERF  building  study  which  concentrates  mostly  on  heat 
movement  within  the  profile  and  total  column  heave.  The  extension  of  the  application  of  the 
model  to  consider  heave  and  potential  for  downslope  movement  is  looked  at  with  the  application 
of  the  model  to  the  CAEN  data  set. 

A  variety  of  means  of  handling  the  data  output  from  the  model  have  been  investigated  and  the 
results  presented  here  are  created  using  cricket  graph  and  UNIRAS. 


1)  Sensitivity  Analysis 

A  sensitivity  analysis  provides  an  opportunity  to  discover  the  most  significant  parameters  in 
determining  model  results  and  also  to  investigate  if  the  model  is  responding  to  changes  in 
parameter  values  in  the  manner  that  one  would  expect. 

To  perform  the  sensitivity  analysis  a  test  profile  was  created  which  was  50cm  deep  and  input 
files  were  used  for  a  silty  loam  soil  type.  A  short  freeze/thaw  cycle  was  used  as  input  to  allow 
for  a  short  simulation  duration  and  permit  many  runs  to  be  made.  With  the  sensitivity  analysis 
made  here  only  one  parameter  was  altered  at  any  one  time.  This  is  very  useful  in  giving  a 
general  idea  of  how  the  model  is  operating  under  different  circumstances  but  means  that 
interrelationships  between  parameter  values  which  inevitably  exist  are  not  considered.  An 
analysis  of  the  combinations  of  timestep  and  cell  size  were  made.  As  cell  size  is  increased 
timestep  can  be  increased.  The  time-step  and  cell  size  limits  are  also  influenced  by  the  hydraulic 
conductivities  set  as  these  have  a  strong  influence  in  controlling  the  mass  transfers  which  will 
cause  instability  if  they  are  too  great  within  any  time-step  length. 

The  functioning  of  the  hydrological  component  of  the  model  is  of  great  importance  as  this  is 
highly  significant  in  determining  water  migration  and  therefore  the  major  influence  on  the 
occurrence  and  extent  of  creep  and  gelifluction  processes.  If  volume  change  in  a  profile  occurred 
only  because  of  the  in  situ  change  of  water  to  ice  then  only  a  small  quantity  of  heave  would 
result.  In  many  cases  the  freezing  of  water  in  the  upper  soil  layers  reduces  the  liquid  water 
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content  of  the  soil  and  creates  a  hydraulic  gradient  sufficient  to  pull  water  up  through  the  profile. 
This  process  continues  until  the  pore  connectivity  reduces  so  that  the  transfer  rates  are  so  slow 
that  movement  is  effectively  halted.  In  certain  silts  ice  contents  can  lead  to  cell  volumes 
doubling. 

As  expected  the  value  of  the  hydraulic  conductivity  was  important.  The  value  for  saturated 
conductivity  is  supplied  to  the  model  and  this  is  then  used  to  create  the  conductivity  with  liquid 
moisture  content  curve  within  the  PERIMOVE  program.  An  increase  in  hydraulic  conductivity 
reduces  the  depth  to  which  the  freezing  front  penetrates  and  increases  the  total  heave  and  the 
excess  ice  contents  in  the  upper  soil  layers.  For  the  test  profile  the  conductivity  range  of  1x10'^ 
ms'^  to  1x10"^  ms'^  produced  the  corresponding  range  in  heave  values  of  33.8mm  to  1.2mm. 
The  slowing  of  frost  penetration  -with  increasing  conductivity  occurs  as  a  result  of  the  increased 
volumes  of  water  in  the  upper  cells  requiring  a  greater  time  period  to  freeze.  Figure  5  illustrates 
the  sensitivity  of  output  to  variations  in  hydraulic  conductivity. 

The  value  of  saturated  moisture  content  was  varied.  The  major  effect  of  this  was  in  the  soil 
heave  values.  As  the  saturated  moisture  content  is  used  to  calculate  soil  porosity  the  effect  of 
increasing  the  saturated  moisture  content  whilst  holding  the  initial  profile  moisture  contents  fixed 
is  to  reduce  heave  as  greater  volumetric  change  in  cell  moisture  content  is  required  before  heave 
can  commence.  As  the  derivation  of  the  conductivity-moisture  curve  utilizes  the  value  of 
saturated  moisture  content  the  conductivities  for  corresponding  moisture  contents  are  reduced 
and  rate  of  water  movement  is  slowed. 

Bulk  densities  (both  frozen  and  dry)  were  altered  although  the  impact  of  altering  these  did  not 
result  in  a  straightforward  relationship.  The  most  significant  impact  was  to  increase  the  depth  of 
frost  penetration  with  increasing  bulk  densities.  Therefore  the  impact  was  mostly  in  the  thermal 
energy  transfer  component. 

The  soil  particle  conductivity  is  very  significant  impact  on  the  thermal  component  of  the  model 
although  not  greatly  influencing  the  hydrological  response.  Increasing  the  thermal  conductivity 
increases  the  rate  and  thereby  the  depth  of  thaw  penetration.  Soil  heave  is  reduced  a  little  which 
probably  reflects  the  slight  reduction  in  water  movement  caused  by  the  profile  freezing  earlier 
and  therefore  water  migration  routes  being  severed  sooner. 

The  sensitivity  analysis  has  identified  the  implications  of  parameter  changes  to  model  operation 
and  has  also  providing  a  guide  to  the  temporal  and  spatial  time-step  limitations  of  the  model. 
This  knowledge  can  be  used  to  aid  the  application  of  the  model  to  laboratory  data  and  test  the 
extend  the  model  validation  to  practical  examples. 


II)  FERF  Building  Simulations 

The  FERF  building  simulations  concentrated  on  the  simulation  of  profile  temperatures  and  total 
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soil  heave  values.  For  the  simulation  a  test  profile  was  created  with  1cm  deep  cells  in  the  upper 
10cm  of  the  profile  and  2cm  deep  cells  in  the  remaining  80cm  to  create  a  total  depth  of  90cm. 
An  example  simulation  is  shown  in  figure  6.  The  parameter  used  to  indicate  the  depth  of  the 
frost/thaw  front  is  the  O^C  isotherm.  With  the  model  this  is  the  temperature  at  which  ice  can 
begin  to  form  but,  as  is  also  true  for  the  observed  situation,  does  not  always  indicate  that  ice 
forms  a  significant  volumetric  portion  of  the  cell. 

The  soil  heave  predicted  by  the  model  is  within  the  range  observed  for  the  study  which  was  from 
2.5cm  to  11cm  and  the  heave  extent  can  be  used  as  an  indicator  of  the  migration  of  water  up 
through  the  profile  to  create  the  observed  ice  contents  and  the  excess  total  water  contents.  The 
simulation  shown  uses  the  data  files  established  for  a  silty  loam  soil  type.  A  certain  amount  of 
model  calibration  was  made  by  adjusting  input  parameter  values  (mostly  thermal  conductivities) 
to  achieve  the  simulation  shown  and  the  basal  temperature  boundary  condition  was  also  used  as 
an  additional  factor  in  the  calibration.  The  depth  of  the  profile  used  in  the  FERF  building  study 
was  Im  with  an  underlying  layer  of  gravelly  sand  of  a  further  2m  and  so  the  lower  zone  of  the 
profile  would  have  acted  as  a  heat  store  to  counter  frost  penetration. 


Ill)  CAEN  Laboratory  Simulations 

The  CAEN  study  provides  additional  sources  of  data  for  model  development  and  validation.  The 
laboratory  investigation  looked  at  the  response  to  freezing  and  thawing  of  four  soil  types.  For 
this  investigation  two  of  those  soil  types  were  selected  to  see  if  the  model  could  successfully 
reproduce  the  differences  between  the  soil  types.  The  physical  properties  of  the  two  soil  types 
selected  are  summarised  in  table  1 .  The  data  was  collected  over  a  suite  of  freezing  and  thawing 
cycles  and  one  of  the  cycles  was  selected  for  the  modelling  program.  The  ambient  temperature 
conditions  and  the  observed  position  of  the  O^C  isotherm  are  shown  on  figure  7.  It  appears  that 
the  sandy  soil  freezes  slower  than  the  silty  soil  which  is  not  what  one  would  expect  however  this 
may  have  been  due  to  the  production  of  large  quantities  of  needle  ice  in  the  upper  few 
centimetres  of  the  profile  which  greatly  slowed  the  initial  frost  penetration.  The  thaw  penetrates 
the  sandy  soil  faster  which  is  as  one  would  expect.  Plots  of  the  observed  and  simulated  depths  of 
the  OX  isotherm,  soil  heave  and  estimates  of  downslope  movement  of  material  are  shown  in 
figure  8  for  the  silt  soil  types.  Figure  9  shows  an  alternative  method  of  displaying  the 
temperature  profiles  over  time  using  a  UNIRAS  program.  Water,  ice  and  other  model  output  can 
also  be  displayed  in  this  way. 

Because  the  CAEN  experiment  used  a  series  of  freeze/thaw  cycles  it  will  be  possible  to  model 
alternative  cycles  with  the  parameter  sets  derived  for  the  cycle  modelled  (cycle  7)  to  give  an 
indication  of  how  much  calibration  of  input  data  is  required  and  if  data  sets  that  reproduce  the 
observed  patterns  over  a  range  of  cycles  can  be  determined.  In  addition  further  work  will  look  at 
the  differences  between  soil  type  responses. 

Heave  values  are  calculated  assuming  soil  expansion  due  to  freezing  occurs  orthogonal  to  the 
slope  and  resettlement  values  by  assuming  gravitational  resettlement  of  the  heaved  material. 


FERF  temperature  input  data 


FERF  simuiated  protile  frost  heave 
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Figure  6:  Example  simulation  of  FERF  building  experiment  using  the  PERIMOVE 
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Table  1;  Summary  of  physical  properties  of  soil  used  in  CAEN  experiment 


CAEN  -  observed  air  temperature,  cycle  7 

Temperature  oC 
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Figure  7:  Observed  ambient  temperature  and  frost  thaw  depths  for  the  CAEN  laboratory 

experiment 
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Observed  and  simulated  results  for  silt  soil  type  CAEN 

Oapth  cm 
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Downslope  movement  at  the  surface  due  to  creep:  0.6  cm 

Downslope  movement  at  the  surface  due  to  gelifluction:  3.6  cm 

Total  downslope  movement  at  surface  for  freeze/thaw  cycle:  4.2  cm 


Figure  8:  Results  from  simulation  of  silt  soil  type  for  the  CAEN  experiment 


Soil  rt*ni|)t‘i;iturt;>  ( 
■■  ABOVE  10 


■  <)■  !  INIR  AS  'iiiiinii  sluiwmi’  profik’  k-ni|n‘r.ituri‘s  over  siiiuilalioii  lor  (.  Al-N  silt  soil 


For  the  CAEN  study  the  slope  angle  was  12°.  It  is  not  really  possible  to  compare  observed  and 
simulated  values  since  the  results  from  the  CAEN  study  were  used  to  determine  relationships 
between  excess  water  contents,  soil  heave  and  movement  due  to  gelifuction.  Further  comparison 
with  other  cycles  and  additional  data  sources  will  help  to  ensure  that  relationships  developed  are 
general  in  nature. 


5  CONCLUSIONS  REACHED 


•  The  PERDVIOVE  model  is  a  complex  model  which  is  capable  of  simulating  the  important 
processes  governing  moisture  and  heat  transfer  through  soils,  soil  heave  and  solifluction. 
The  model  operation  and  responses  to  changes  in  parameters  and  input  data  are  as 
expected. 

•  very  good  simulations  of  observed  data  have  been  achieved  with  the  model  for  a  range  of 
experimental  conditions. 

•  Data  collection  in  soil  freezing  and  thawing  environments  is  difficult  but  new  techniques 
are  being  developed  and  many  of  these  concentrate  on  methods  for  measuring  liquid 
moisture  contents  within  the  soil  environment.  The  collection  of  further  'internal'  data 
will  place  further  constraints  on  model  validation  and  ensure  that  the  model  is  behaving 
in  the  desired  manner.  Routes  for  further  interaction  with  experimental  projects  have 
been  identified. 

•  The  PERIMOVE  model  can  also  be  viewed  as  a  tool  to  investigate  the  influence  of 
changing  the  controlling  factors  governing  freezing,  thawing  and  mass  movement  to  asses 
the  influence  of  these. 


6  FUTURE  LINES  OF  RESEARCH  ARISING  FROM  THE  PROJECT 

Two  main  lines  of  further  research  can  be  identified  which  would  enhance  the  capability  and 
applicability  of  the  model.  They  divide  into  further  improvement  to  the  structure  of  the  model 
and  extensions  to  the  testing  and  application  of  the  model  in  association  with  experimental 
results. 

Firstly  to  be  considered  is  the  extension  to  a  fully  2-dimensionaI  scheme  incorporating  a 
constantly  deforming  grid  element  structure.  These  are  highly  complex  to  program  as  model 
equations  must  then  incorporate  a  term  which  includes  the  velocity  due  to  the  movement  of  the 
grid  as  well  as  the  velocity  due  to  the  heat  or  water  exchanges.  The  structure  of  such  schemes 
has  been  investigated  and  a  connection  made  with  Mary  Albert  at  CRREL  (see  for  example 
Albert,  1984  and  Albert  and  O'Neill  1986).  Figure  10  shows  diagrams  of  a  deforming  mesh  from 
Albert  (1984)  where  the  movement  of  the  phase  boundary  was  modelled  for  a  drum  filled  with 
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saturated  sand  in  which  freezing  was  instigated  by  the  passage  of  cold  ethylene  glycol  through  a 
vertical  copper  pipe  situated  off-centre  within  the  drum.  It  may  be  possible  that  with  the 
continued  development  of  finite  element  programming  and  automatic  mesh  generation 
techniques  it  wiU  eventually  be  possible  to  obtain  standard  deforming  grid  packages  to  which  the 
model  can  be  linked.  Alternatively,  the  connections  established  at  CRREL  may  provide  the  basis 
for  the  development  of  a  scheme  to  directly  link  in  with  the  current  model. 

Secondly,  as  mentioned  in  the  previous  section  is  the  desire  to  enhance  the  model  by  further 
interaction  with  observed  data  sources.  Harris  et  al  (1993)  are  continuing  the  project  based  at 
CAEN  with  a  further  series  of  freeze/thaw  experiments  planned  in  1994.  Additional  data  on  soil 
moisture  status  during  freezing  will  be  collected  then  and  the  aim  is  to  collect  a  more 
comprehensive  data  set  particularly  with  regard  to  internal  profile  conditions  during  freezing  and 
thawing.  They  aim  to  concentrate  on  two  different  soil  types  to  highlight  the  differences  between 
coarse  grained  sandy  soils  and  finer  grained  silty  soils. 

In  addition  to  the  use  of  the  model  to  simulate  specific  experimental  conditions  it  should  be 
possible  to  extend  the  use  of  the  model  as  a  tool  in  investigating  controlling  factors  on  slope 
movement  and  development  in  cold  climates.  This  could  be  to  investigate  movements  in  current 
cold  regions  or  to  try  to  unravel  climatic  conditions  which  could  have  led  to  the  formation  of 
certain  geomorphological  features. 
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-  frozen  soil 

-  unfrozen  soil 


CONDITION  (1) 

critical  depth  of 
frozen  soil  which  will 
support  chosen  vehicle 


CONDITION  (2) 

critical  depth  of  unfrozen 
(’slippery’)  layer  where  traction 
is  too  low.  If  tires  reach 
frozen  layer  then  zone  may  be 
trafficable. 


CONDITION  (3) 

soil  moisture  status  is  critical 
for  determining  trafficability. 

The  frozen  layer  is  too  far 
down  the  profile  to  give  support 
and  acts  to  impede  drainage. 


Figure  1 1;  Summary  diagram  of  critical  soil  trafficability  conditions  in  freeze/thaw 

environments. 


APPENDIX  1 


Incut  files  required  to  run  PERIMOVE  and  parameter  definitions. 


GUIDE  TO  THE  INPUT  FILES  REQUIRED  TO  RUN  PERIMOVE 


The  input  files  are  listed  in  alphabetical  order.  All  inputs  are  in  free  format  unless 
otherwise  stated. 


INDEX: 

Ambientl.data 

Ambt.data 

Cvin.data 

Hydro.data 

InitialH.data 

InitialT.data 

Kcurve.data 

Out.data 

SM.data 

Soli,  data 

Struc.data 

TW.data 

Temporal. data 

Ttsec.data 


initial  ambient  temperature  and  number  of  temperature 
changes  in  simulation. 

array  of  ambient  temperatures  for  simulation. 

soil  freezing  parameter  data. 

soil  hydrological  property  data. 

soil  hydrological  condition  at  start  of  simulation. 

soil  temperature  condition  at  start  of  simulation. 

ice  charateristic  curve  data. 

output  option  controls. 

data  for  suction/moisture  curve. 

data  for  solifluction  potential. 

cell  configuration  and  definition  of  computational  points, 
data  for  temperature/water  curve, 
controls  for  length  of  simulation, 
time  of  each  temperature  change. 


Ambientl.data  (unit 49) 


LI 


NOTC,  STAMBT,  BTEMP 

NOTC  number  of  temperature  changes  in  mn. 
STAMBT  initial  ambient  temperature  (K). 

BTEMP  basal  boundary  temperature  (K). 


Ambt.data  (unit  48) 

LI  AMBT(NOTC) 

AMBT  ambient  temperature  (K). 


L1  is  repeated  NOTC  times. 


Cvin.data  (unit  59) 


LI  UDENS,  FDENS,  SHCSOIL,  CW,  Cl 

UDENS  unfrozen  soil  bulk  density  g  m'^. 

FDENS  frozen  dry  bulk  density  g  m‘3. 

SHCSOIL  specific  heat  capacity  of  soil  J  g"^  K‘^ . 

CW  specific  heat  capacity  of  soil  water  J  g"^  K‘^ . 

Cl  specific  heat  capacity  of  ice  J  g"^  K'^ . 

Hydro.data  (unit  57) 

L1  STCON1.STCON2 

STCON1  saturated  soil  hydraulic  conductivity  for  soil  layer  1  m  s'^ 

STCON2  saturated  soil  hydraulic  conductivity  for  soil  layer  2  m  s’^ . 

LJ2  ASR1,ASR2 

ASR1  saturated  soil  moisture  content  for  soil  layer  1  m'^. 

ASR2  saturated  soil  moisture  content  for  soil  layer  2  m^  m'^. 

L3  (HYDG(I),  l=1,NC) 

HYDG  hydraulic  gradient  under  saturated  flow  -  based  on  slope  angle. 
NC  number  of  columns  in  simulation. 

L4  (ORESt.M.NC) 

ORES1  residual  soil  moisture  content  for  layer  1 . 

L5  (ORES2,  l=1,NC) 

ORES2  residual  soil  moisture  content  for  layer  2. 

InitialH.data  (unit  54) 

L1  (ATHESP(I,J).J=1,M) 

ATHESP  initial  soil  moisture  contents  m^  m'^. 

M  number  of  ceils  in  column. 

•  L1  is  repeated  NC  times. 


f 

iii 

L2 

(VATHET(I.J).  J=1,M) 

VATHET  ice  content  or  potential  volume  of  moisture  available  for  release 
on  thaw  m'^. 

InitialT.data 

(unit  53) 

L1 

(sltp(I.J).  J=1.M) 

sitp  initial  soil  temperature  (K)  for  each  cell. 

M  number  of  cells  in  column. 

*  L1  is  repeated  NC  times  (1=1  ,NC). 

Kcurve.data 

(unit  58) 

L1 

NKI 

NKI 

number  of  values  on  thermal  conductivity  of  ice  with  temperature 

curve  (Ktl/T). 

L2 

Kl 

(KI(I).I=1,NKI) 

thermal  conductivity  of  ice  on  Ktl/T  curve  W  m'^  K"^ . 

L3 

TKI 

(TKI(I).  1=1. NKI) 

temperature  values  on  Ktl/T  curve  K. 

L4 

NKO 

NKO 

number  of  values  on  thermal  conductivty  of  water  with  temperature 
curve  (Kt0/T). 

L5 

KO 

(KO(l),  1=1, NKO) 

thermal  conductivity  of  water  on  Kt0/T  curve  W  m'^  K'^ . 

L6 

TKO 

(TKO(I),  1=1, NKO) 

temperature  values  on  Kt0/T  curve  °K. 

L7 

TDEP 

temperature  at  which  ice  starts  to  melt 

L8 

(KS(I,J),J=1,NNL(I)) 

KS  intrinsic  thermal  conductivity  of  soil  particles  W  m'^  K'^ . 

IV 

NNL  number  of  cells  in  column 

*  L8  repeated  for  each  column  (1=1  ,NC). 

Out.data  (unit  61) 

LI  outmin,  nops 

outsec  iteration  period  for  output  of  data  (minutes), 
nops  number  of  output  data  options. 

L2  (option(l),  1=1, nops) 

option  output  data  option  codes,  (see  key  below) 

L3  nslices 

nslices  number  of  columns  for  which  graphics  data  are  to  be  written  out. 

L4  (slice(K),  K=1, nslices) 

slice  column  for  which  graphics  data  is  to  be  written  out. 

L5  NOCOP  number  of  cells  for  which  output  is  required 

L6  (CELLOP(l),l=1, NOCOP) 

cell  numbers  for  which  output  is  required 
*L5  and  L6  are  repeated  for  each  column  (NC  times) 

L7  outstart  minute  from  which  output  is  produced 


Key  to  PERIMO  output  options. 

Numerical  data  files  with  captions. 


OPTION 

FILE  NAME 

1 

sitp.data 

2 

ice.data 

3 

water.data 

4 

nhflux.data 

5 

hflux.data 

6 

heatcap.data 

7 

thermcon.data 

8 

hydcon.data 

9 

suction.data 

10 

fcum.data 

11 

fmax.data 

12 

fmin.data 

13 

fs.data 

FORTRAN 

UNIT  NO. 

DATA 

81 

slip 

82 

VATHET 

83 

ATHESP 

84 

NHFLUX 

85 

HFLUX 

86 

VHC 

87 

APKT 

88 

KU 

89 

set 

90 

FCUM 

91 

FMAX 

92 

FMIN 

93 

FS 

14 

athi  data 

94 

ATHI 

15 

vathi  data 

S5 

VATHI 

16 

atempdata 

96 

ATEMP 

» 

btf!ux.data 

97 

BTFLUX 

40 

move.data 

40 

GDISP+HDISP 

Graphical  data  files 

OPTION 

FILE  NAME 

FORTRAN 
UNIT  NO. 

DATA 

20 

Gsitp.data 

70 

sitp 

21 

Gice.data 

72 

VATHET 

22 

Gwater.data 

73 

ATHESP 

23 

Gsuct ion  .data 

73 

set 

24 

Ghycoixl.data 

74 

KP 

25 

gtemp.data 

75 

ATEMP 

26 

Gkt.data 

76 

KT 

27 

Gapkt.data 

77 

APKT 

28 

Gfs.data 

78 

FS 

29 

Gnhflux.data 

79 

NHFLUX 

30 

Ghflux.data 

80 

HFLUX 

33 

Gexp.data 

33 

COLHT,  COLHT-DEPTHP 

34 

Gcelexp.data 

34 

EXP 

35 

fort.300-400 

(isotherm  position)  ISO, KKK 

e.g.  fort.302  contains  ISO  ar>d  KKK  data  for  column  2 

36 

Gsoilp.data 

36 

SOILP 

37 

Gexwat.data 

37 

EXWAT 

38 

Gdisplace.data 

38 

HDISP.  GDISP 

SM.data  (unit  56) 


NQ1 

NQ2 

KQ 


NQ1.  NQ2,  KQ 

number  of  observations  on  suction/moisture  curve  for 
hydrological  layer  1 . 

number  of  observations  on  suction/moisture  curve  for 
hydrological  layer  2. 

number  of  observations  on  hydraulic  conductivity/moisture  curve. 


L2 


X1 


(X1(I),I=1.NQ1) 

moisture  values  on  suction/moisture  curve  for  layer  1  m^  m'^. 


L3 


Y1 


(Y1(I).I=1.NQ1) 

suction  values  on  suction/moisture  cun/e  for  layer  1  m. 


lines  L2  and  L3  are  repeated  for  layer  2  for  X2  and  Y2. 


Soli.data  (unit  60) 


L1  GAMMAS,  GAMMAU 

GAMMAS  saturated  unit  weight  of  soil  g  m'^. 
GAMMAU  unsaturated  unit  weight  of  soil  g  m"^. 

L2  (BETA(I),  l=1.NC) 

BETA  slope  angle  in  degrees 

L3  (C{I,J),  J=1.M) 

C  cohesion  KPa  for  each  cell  (J=1  ,M). 

*  L3  is  repeated  for  each  column. 

L4  (PHKI.J),  J=1,M) 

PHI  internal  friction  angle  in  radians. 

*  L4  is  repeated  for  each  column. 

L5  GELRAT  heave:gelif!uction  ratio  for  soil  type 

Structural. data  (unit  51) 

L1  NC,  NOC,  LENGTH,  BDT 

NC  number  of  columns  in  simulation 

NOC  number  of  cells  in  first  column  of  thp  p'-ofile 

LENGTH  cell  length  (m). 

BDT  breadth  of  cell  (m) 

L2  (NNL(I),  l=1,NC) 

NNL  number  of  cells  in  each  column  of  the  profile. 

L3  (CDEPTH(I,J),  dcp{l,J),  1=1  ,NC,  J=1  ,NNL(I)) 

CDEPTH  cell  depth  (m) 

dcp  depth  to  computation  point  in  cell(m) 

L4  (DEPTHP(I),  1=1  ,NC) 

DEPTHP  depth  of  profile  in  each  column  (m). 


L5 


(HT(I,J),  J=1,M) 


HT 


height  of  each  computation  point  above  the  profile  base. 

*  L3  is  repeated  for  each  column 
L6  (LC1(1),  l=1.NC) 

LC1  position  (j)  of  last  cell  in  hydrological  layer  1 . 

L7  (FC2(I).  I=1.NC) 

FC2  position  (j)  of  first  cell  in  hydrological  layer  2. 

L8  (TTCOM{I.J),  J=1.M) 

depth  of  cell  computation  point  from  surface  (m). 

'  line  L8  is  repeated  for  each  column. 

TW.data  (unit  55) 

L1  TF 

TF  number  of  observations  on  characteristic  freezing  curve 

L2  (TEMCV(I),  l=1,TF) 

TEMCV  temperature  values  on  characteristic  freezing  curve,  K  (highest 
first). 

L3  (WATECV(I),I=1,TF) 

WATECV  moisture  values  on  characteristic  freezing  curve  (m^  m'^). 

Temporal.data  (unit=50) 

L1  IP,  NOSECR 

IP  iteration  period  (seconds). 

NOMINR  number  of  minutes  in  simulation. 

Ttsec.data  (unit=52) 

L1  TTSEC(JJJ) 

TTSEC  time  from  start  of  simulation  to  temperature  change  input  to 
PERIMO  model  (minutes). 

number  of  temperature  changes  during  simulation  (NOTC). 


JJJ 


DEFINITIONS  OF  VARIABLES  FOR  PERIMOVE 


AMBT . ambient  temperature  K. 

APKT . apparent  thermal  conductivity  W  m’''  K'^  . 

ASR1 . saturated  moisture  content  for  soil  type  1  m^  m‘^. 

ASR2 . saturated  moisture  content  for  soil  type  2  m^  m'^. 

ATEMP . Ambient  temperature  (K). 

ATH1 . volume  of  water  in  cell. 

ATHESP . soil  water  content  m^  m'^. 

ATHESP1 . soil  water  volume  m^. 

ATHET . volumetric  water  content  derived  from  T/0  curve. 

ATHET 1 . water  content  derived  from  T/0  curve,  m^. 

AVK . averge  hydraulic  conductivity  m  s‘^ . 

AVKT . Average  thermal  conductivity  across  cell  boundary. 

BDT . cell  breadth  m. 

BETA . slope  angle. 

BFLUX . volume  of  flux  through  base  of  profile. 

BISO . depth  of  basal  isotherm. 

BWFLUX . basal  flux. 

C . cohesion. 

CDEPTH . Depth  (size)  of  ceils  (m). 

CELVOL . cell  volume. 

CH . Temperature  change. 

CHPERC . Temperature  change  expected  given  net  input  or  loss  of  energy  on  phase 

change.  Calculated  per  cell. 

COLHT  height  of  each  soil  column. 

cmdt . Thermal  component  of  energy  budget  for  dT. 

CVGR . curve  gradient  on  T/0  curve. 

dcp . depth  to  computation  point  in  cell. 

deltaT . Change  in  temperature  during  iteration  period. 

deltaw . Change  in  moisture  during  iteration  period. 

DEPTHP . depth  of  profile  in  each  column. 

DISC . depth  from  surface  at  which  273.15K  isotherm  occurs. 

EXP . expansion  of  cell  in  vertical  from  original  volume. 

EXSAT . change  in  moisture  content  due  to  saturated  flow. 

EXWAT . excess  water  content  expelled  from  cell  on  thawing. 

FC2 . position  0)  of  first  cell  in  layer  2. 

FCUM . .cumul?*ive  factor  of  safety  per  column. 

FDENS . frozen  dry  bulk  density  (g  m'^). 

FDMAX . depth  of  max.  FS. 

FDMIN . depth  of  min.  FS. 

FLUX . soil  water  flux  between  ceils,  volume. 

FMAX . maximum  factor  of  safety  per  column. 

EMIN . minimum  factor  of  safety  per  column. 

FS . factor  of  safety. 

G1,  G2 . average  incremental  gradients  on  suction/moisture  curve,layers  1  and  2. 

GAMMAS . saturated  unit  weight  of  soil. 

GAMMAU . unsaturated  weight  of  soil. 

GZ1 ,  GZ2 . average  incremental  gradients  on  K/0  curve. 

HFLUX . Heat  flux  across  cell  boundary. 

HPOT . hydraulic  potential. 

HT . height  of  each  computation  point  above  profile  base. 

HYDG . hydraulic  gradient  sat.  flow,  ie.  slope  angle. 

IKG . incremental  gradient  on  Ktl/T  curve. 

IP . iteration  period  secs. 

•SO . depth  from  previous  compn.  point  at  which  273.15K  isotherm  occurs. 


ITT . count  for  iterations. 

Kl . thermal  conductivity  of  ice,  Ktl/T  curve. 

KIK . thermal  conductivity  of  ice  at  given  temp. 

KO . thermal  conductivity  of  water,  KtO/T  curve. 

KOK . .thermal  conductivity  of  water  at  given  temp. 

KP . cell  hydraulic  conductivity,  m  s'"’ . 

KQ . number  of  observations  on  K/0  curve. 

KS . intrinsic  thermal  conductivity  of  soil  solids. 

LA . Latent  heat  of  fusion. 

LAO . hydrological  component  of  energy  release  for  dT. 

LC1 . position  (j)  of  last  cell  in  layer  1 . 

LENGTH . cell  length,  m. 

MORES . Mean  residual  moisture  contents. 

NO . no.  of  columns. 

NEWVOL . updated  cell  volume  followind  expansion  or  resettlement. 

NFLUX . net  water  flux  of  cell. 

NHFLUX . net  heat  flux  W  m‘^  s‘^. 

NKI . number  of  values  on  Ktl/T  curve. 

NKO . number  of  values  on  KtO/T  curve. 

NNL . number  of  cells  in  each  column  of  profile. 

NNLMAX . maximum  number  of  cells  in  any  of  the  columns. 

NOC . number  of  cells  in  first  column  of  the  profile. 

NOIT . number  of  iterations  in  simulation. 

NOMINR . number  of  minutes  in  run. 

NOSECR . number  of  seconds  in  run. 

NOTC . .number  of  temperature  changes. 

NQ1,  NQ2 . number  of  observations  on  suction/moisture  curve,  layers  1  and  2. 

NSEIT . number  of  seconds  passed  by  end  of  iteration. 

nslices . number  of  columns  for  which  graphics  data  are  written. 

OKG . incremental  gradient  of  KtO/T  curve. 

option . output  data  options. 

ORESt,  ORES2  .  .  .  residual  soil  moisture  content,  layer  1  &  2. 

ORIVOL . cell  volume  at  start  of  simulation. 

OUTMIN . iteration  period  for  output  data,  minutes. 

OUTSEC . iteration  period  for  output,  minutes. 

OUTSTART . control  for  start  time  for  output,  minutes. 

PHI . internal  friction  angle  of  soil. 

set . soil  suction  in  each  cell. 

SEC . difference  between  nseit  and  ttsec  input. 

SFI . shape  factor  ice. 

SFO . shape  factor  water. 

SFS . shape  factor  solids. 

SHCSOIL . specific  heat  capacity  of  soil  (J  m'3  K'^). 

NC . number  of  columns. 

slice . column  for  which  graphics  data  is  to  be  output. 

sitp . soil  temp.  (K). 

sltpt2 . soil  temp  calculated  by  heatrans.f  on  basis  of  changing  ambient  temp. 

SOILP . percentage  volume  of  swi  particles  per  cell. 

SOILV . volume  of  soil  solids,  m'^. 

STAMBT . starting  ambient  temperature. 

STCH . ambient  temperature. 

STCON1  STCON2  .  saturated  hydraulic  conductivity  m  s'^ . 

TEMCV . temperature  values  on  T/0  curve. 

TEMCV . Temperature  values  on  T/0  curve  (K). 

TF . number  of  observations  on  T/0  curve. 

tgrad . temperature  gradient  across  frozen  fringe. 

TKI . temperature  values  om  Ktl/T  curve. 

TKO . temperature  values  on  KtO/T  curve. 


TQ . /lumber  of  observations  on  Ts/0  curve. 

TTCOM . .depth  of  cell  computation  point  from  surface. 

TTMIN . time(minutes)  from  start  of  run  to  temperature  change. 

TTSEC . time  from  start  of  run  for  temperature  change. 

UDENS . unfrozen  dry  bulk  density  (g  m-3). 

VADFAC . volume  adjustment  factor  for  cell  expansion/contraction. 

VATH1 . volume  of  ice  in  cell. 

VATHET . ice  content  m3  m-3. 

VOLCH . ratio  of  volume  to  original  cell  volume. 

VHC . .volumetric  heat  capacity  of  cell  y  m"^  K'^ ). 

WATECV . Moisture  values  on  T/0  curve  m^  m'^. 

XI ,  X2 . moisture  values  on  suction  moisture  curve,  layers  1  &  2. 

Y1,  Y2 . suction  values  on  suction  moisture  curve,  layers  1  &  2. 

Z1 ,  Z2 . conductivity  values  on  conductivity  moisture  curve,  layers  1  &  2. 


APPENDIX  2 


PERIMOVE  FORTRAN  CODE 


PERIMOVE.f 


Cv.f . j(ii 

freeze.! . xiv 

heatrans.f . xviii 

isotherm.! . xxi 

kta.! . xxiv 

ku.! . xxvi 

masstrans.! . xxvii 

output.! . xxxi 

soli.! . xW 

thaw.! . xliii 


c  Program  PERIMOVE.f 

c 

^  «•*************•**%*«*»*«***«••**•«•»*«*********«**•**«***••*«*«* 

c  PERIMOVE.F 

^  .......................................................... 

c  Program  to  predict  the  volume  of  ice  present  in  a  soil  profile  during 

c  freezing  and  thawing.  The  program  also  predicts  heave  and  resettlement 

c  of  the  soil  profile  and  volumetric  water  excesses  following  thawing, 
c 

c  The  program  calculates  the  soil  temperature  at  depth,  through 
c  time,  modifies  the  temperature  for  transfer  of  latent  heat 

c  and  updates  soil  moisture  content  accordingly. 


Q  ******************************************************************** 

c  program  PERI 

program  perimove 

^  ******************************************************************** 

c  DECLARATION  OF  VARIABLES 

^  *'**«««•*«**************«**«***«**«***«************«******«********** 

DOUBLE  PRECISION  TTCOIVI(5,60),  HT(5,60),COLHT(5) 

DOUBLE  PRECISION  LENGTH,  CDEPTH(5,60),  BDT,CELVOL(5,60) 
DOUBLE  PRECISION  dcp(5,60),ORIVOL(5,60),VADFAC(5,60) 
DOUBLE  PRECISION  DEPTHP(5),SOILP(5,60),EXP(5,60),CDEP(5,60) 
DOUBLE  PRECISION  ATHET(5,60),  VATHET(5,60),ATHESP(5,60) 
DOUBLE  PRECISION  ORES1(50),  ORES2(50),  XI  (50),  X2(50) 
DOUBLE  PRECISION  Y1(50),  Y2(50)  ,STCON1,  STCON2 
DOUBLE  PRECISION  ASR1,  ASR2,  sltp(5,60),  sltpt2(5,60) 

DOUBLE  PRECISION  KP(5,60),sct(5,60),EXWAT(5,60) 

DOUBLE  PRECISION  LA,  VHC(5,60),  G1(50),  G2(50) 

DOUBLE  PRECISION  Z1(50),  Z2(50),  GZ1(50),  GZ2(50) 

DOUBLE  PRECISION  FLUX(5,60),  NFLUX(5,60),  HPOT{5,60) 

DOUBLE  PRECISION  BWFLUX(50),  STATH(5,60),  STVATH(5,60) 
DOUBLE  PRECISION  AVK(5,60),  HYDG(50) 

DOUBLE  PRECISION  TEMCV(50),  WATECV(50),  CVGR(50),TDEP 
DOUBLE  PRECISION  APKT(5,60),  KT(5,60) 

DOUBLE  PRECISION  HFLUX(5,60),  NHFLUX(5,60) 

DOUBLE  PRECISION  AMBT(525600),  ATEMP(1 468800) 

DOUBLE  PRECISION  STAMBT,BTEMP 
DOUBLE  PRECISION  FS(5,60),  FMIN(50),  FMAX(50) 

DOUBLE  PRECISION  FDMAX(50),  FDMIN(50),  FCUM(50) 

DOUBLE  PRECISION  EXSAT(50).1SO(5,60),SOIL1  (5,60) 

DOUBLE  PRECISION  SHATHET(5,60),  SHVATHET(5,60) 

DOUBLE  PRECISION  HTHETAC(50),KI(20),TKI(20),KS(5,60) 

DOUBLE  PRECISION  KO(20),TKO(20),OKG(20),IKG(20) 

DOUBLE  PRECISION  C(5,60),  BETA(5),MAXCOLHT(5),HDISPI5) 
DOUBLE  PRECISION  GAMMAS,  GAMMAU,GELRAT,GDISP(5) 
DOUBLE  PRECISION  PHI(5,60) 


c  Integers 

INTEGER  NC,  NNL(5),  NOC 
INTEGER  NNLMAX,  col(50),  ncol 
INTEGER  SEC(525600) 
INTEGER  outsec 
real  outmin 

INTEGER  NQ1,NQ2,KQ,TF 
INTEGER  LC1(50),  FC2(50) 
INTEGER  option(50) 

INTEGER  IP 


i  i 


INTEGER  NOTC 

INTEGER  NSEIT(1 468800),  IH,  NOIT 

INTEGER  NOSECR,NOMINR,TTSEC(525600),nMIN(525600) 

INTEGER  slice(50),  nslices 

INTEGER  isocol(50),  isoncol 

INTEGER  NKO,  NKI.KKK 

INTEGER  NOCOPV,NNLMA,CELL 

INTEGER  NOCOP(30),CELLOP(5,60) 

INTEGER  OUTSTART.OUTTIM 

COMMON  DUMMY(5,60) 
c 

^  ******************************************************************* 

c  OPEN  STATEMENTS 

-  ******************************************************************* 

open(unit=20,  file='DATAIN/Mb.clata',  status='old’) 
rewind  20 

open(unrt=21,  file='DATAIN/lso.dafa‘,  status='old') 
rewind  21 

c  open(unil=32,  file='DATAIN/massbal.data',status=’unknown') 

c  rewind  32 

open(unit=48,  fiie='DATAIN/Ambt.data',  status='old') 

open(unif=49,file='DATAIN/Ambient1.data',sfatus='old') 
rewind  49 

open(unit=50,  fiie='DATAIN/Temporal.data‘,status='old') 
rewind  50 

open(unit=51,  file='DATAIN/Struc.data',  status=‘old') 
rewind  51 

open(unit=52,file=’DATAINn'tsec.data',status='old') 

open(unit=53,file='DATAIN/lnitialT.data',  status='old') 
rewind  53 

open(unit=54,file='DATAIN/lnitialH.data',status='old’) 
rewind  54 

open(unit=55,file='DATAIN/TW.data'.status='old‘) 
rewind  55 

open(unit=56,file='DATAIN/SM.data',status='old’) 
rewind  56 

open(unit=57,file='DATAIN/Hydro.data',status=‘old') 
rewind  57 

open(unit=58 ,  file='DATAIN/Kcunre.data‘,  status='old') 
rewind  58 


open(unit=59,  file='DATAIN/Cv.data',  status='unknown') 
rewind  59 


open(unit=60,  file='DATAlN/Soli.data',status='old') 
rewind  60 

open(unit=61 ,  file='DATAlN/Out.data',  status='old') 
rewind  61 

OUTTIM=1 

***«*********«*«****t*«**********t •«****«•**•***«*••*•«•****«******** 

TEMPORALdata  unit  50 


Read  iteration  time  and  no.  of  minutes  in  mn. 
read(50,*)lP,  NOMINR 

convert  no.  of  minutes  in  run  to  no  of  seconds 
NOSECR=NOMINR*60 

*************«***************«***«***«***«**«**********«**««**t*****t 

STRUCTURAL.dafa  unit  51 

********************************************************************* 

Space  dimensions. 

READ(51,‘)NC,  NOC,  LENGTH,  BDT 

Number  of  cells  in  each  column 
READ(51,*)(NNL(I),I=1,NC) 

calculate  volume  of  each  cell  and  store  in  ORIVOL  for  volume 
change  calculations.  Initialise  volume  adjustment  factors. 

DO  16 1=1, NC 
D017J=1,NNL(1) 

READ(51,*)  CDEPTH(l,J),dcp{(,J) 

CELVOL(l,J)=LENGTH'BDT*CDEPTH(l,J) 

ORIVOL(I.J)=CELVOL(l,J) 

VADFAC(I,J)=1.0 

EXP(l,J)=0.0 

CDEP(I,J)=CDEPTH(I,J) 

CONTINUE 

CONTINUE 


NNLMAX=0 
D01  1=1, NC 

IF(NNL(I).GT.NNLMAX)THEN 

NNLMAX=NNL(I) 

ENDIF 

CONTINUE 

Depth  of  profile  in  each  column. 

READ(51  ,*)(DEPTHP(I),I=1  .NC) 

Height  of  computation  point  above  profile  base. 

D02I=1,NC 

M=NNL(I) 

READ(51,*)(HT(I,J),J=1,M) 

CONTINUE 

Define  LC1  and  FC2. 

READ(51,*)(LC1(I),1=1,NC) 

READ(51,*){FC2(1),I=1,NC) 


IV 


c 

c  Depth  of  each  computation  point  from  surface. 

DO  3  1=1, NC 
M=NNL(1) 

READ(51,*)(TTCOM(l,J).J=1,M) 

3  CONTINUE 
do  3012  1=1, NC 
write(39,*)depthp(i) 

MAXCOLHT(l)=DEPTHP(l) 
do  3013J=1,M 

WRITE(39,‘)l,J,HT(l,J),TTCOM(U).CELVOL(l,J) 

3013  CONTINUE 
3012  CONTINUE 

^  vacatvtc********************************************************** 

c  INITIALT.data 

^  *•**«**• •••**************««**•«•**•****«*•• • > •  •  ••*•*■•«*********•* 

V. 

c  Starting  soil  temperatures 

DO  4  1=1,  NC 
M=NNL(I) 

READ(53,*)(sltp(l,J),J=1,M) 

4  CONTINUE 

^  *****t****t*****************t***»#»***»**»********************»*** 

c  INITIALH.data  unit  54 

^  **************»****»»*******t***t*«*»***»**t*»t»**t**»*****»*#*** 

c  Starting  moisture  conditions 

DO  5  1=1, NC 
M=NNL(I) 

read(54,*)(ATHESP(l,J),J=1,M) 

5  continue 

c  Starting  ice  content 

DO  60 1=1, NC 
M=NNL(I) 

read(54,*)(VATHET(l,J),J=1,M) 

60  continue 

^  ♦***********•♦♦♦♦**♦♦**♦•♦**••********♦*♦*************♦*********** 

c  TW.data  unit  55 

^  **•**■•  ************************************************************ 

c  Nun.ber  of  values  on  soil  water  characteristic  curve  for  temps 
c  below  273. 1 5K.  (Temperature  and  moisture  values) 

c  Ambient  thermal  conductivity 

READ(55,*)TF 

c  Values  on  soil  water  characteristic  curve. 

READ(55,*)(TEMCV(I),I=1,TF) 
READ(55,*)(WATECV(I),I=1,TF) 

^  t**************************************************************** 

c  Ktcurve.data  unit 

Q  **«*******•***«******•******«••«*•******•********•*•******•*«**** 

c  Read  no.  of  values  on  Kl/T  curve 

read(58,*)NKI 

c  Read  Kl  and  T  values  on  Kl/T  curve 

read(58,*)(KI(l),l=1.NKI) 
read(58,*)(TKI(l),l=1,NKI) 


V 


Read  no.  of  values  on  KO/T  curve 
read(58,*)NKO 

Read  KO  and  T  values  on  KO/T  curve 

read(58.*)(KO(l).!=1,NKO) 

read(58,*)(TKO(l),l=1,NKO) 

Read  thaw  initiation  temperature  TDEP 
read(58,*)TDEP 

Read  values  for  intrinsic  thermal  conductivity  of  soil  particles 
do  25 1=1, NC 

read(58,*)(KS(l,J),J=1,NNL(l)) 

continue 


Establish  thermal  conductivity  of  water 
Incremental  gradients  on  K/0  curve 
do  26  l=1,NKO-1 

OKG(l)=(KO(l+1)-KO(l))/(TKO(l+1)-TKO(l)) 

continue 

Incremental  gradients  on  K/l  curve 
do  27  l=1,NKI-1 

IKG(I)=(KI(I+1)-KI(I))/(TK1(I+1)-TKI(I)) 

continue 


•*•••*•**#**••**•**••***•««««***«*«*•«*«#*****«*««****«*********** 

SM.dala  unit  56 

*********»****»**»***tt***»*iHt*****t*»***t**********»t**«t******** 

Number  of  observations  on  suction  moisture  curve,  Ku/moisture 
curve. 

READ(56,*)NQ1,  NQ2,  KQ 


Values  of  moisture  and  suction  on  suction  moisture  curve. 
READ(56,*)(X1(I).I=1,NQ1) 

READ(56,*)(Y1(I),I=1,NQ1) 

READ(56,*)(X2(I),I=1,NQ2) 

READ(56,*)(Y2(I),I=1,NQ2) 

t*****t**t*****i^****t*tt*t****t***t^********t***i^**^t*****tt**t******t 

HYDRO.data  unit  57 

***********•«***«***«««*•••««**«**•****•**«*********'•**««*******'***** 

Saturated  hydraulic  conductivity 
READ(57,*)STCON1,  STCON2 

Saturated  moisture  content  for  layers  1  and  2,  and  hydraulic  gradient 
(slope  angle) 

READ(57,*)ASR1,ASR2 

READ(57,*)(HYDG(I),I=1,NC) 

Initialise  starting  volumetric  contents  of  soil  solids 
DO  591 1=1, NC 


M=NNL(I) 

DO  592  J=1,M 

if(J.LE.LC1(l))SOILP(l,J)=1.0-ASR1 

if(J.GE.FC2(l))SOILPC,J)=1.0-ASR2 

SOIL1(l,J)=SOILP(l,J)'CELVOL(l,J) 

592  continue 

591  continue 

c 

c  Residual  soil  moisture  content  for  layers  1  and  2. 

READ(57,*)(ORES1(l),l=1,NC) 
READ(57,*)(ORES2(l),l=1,NC) 

^  •«***«■««*******•«*«#***«*•****•«**«**«**•«***•••***»*«********«» 

c  READ  IN  STABILITY  DATA 

^  *******************************  *«**»**«*•**«**•**•*****•••«**«••* 

READ(60,*)GAMMAS,  GAMMAU 

READ(60,*)(BETA(I),I=1,NC) 

DO  40  1=1, NC 
M=NNL(I) 

READ(60.*)(C(I,J),J=1,M) 

40  CONTINUE 

DO  41  1=1, NC 
M=NNL(I) 

READ(60,*)(PHI(I,J),J=1,M) 

41  CONTINUE 

DO  42  1=1, NC 
M=NNL(I) 

READ(60,*)GELRAT 

42  CONTINUE 

^  *««**««**«***«*****«*****«*«*«******«**«**«********«t*ft ********** 

c  OUT.data  unit  61 

^  ***************************************************************** 

c  Iteration  period  of  output,  number  of  data  sets. 

read(61,’)outmin,nops 
c  convert  iteration  period  to  seconds 

outsec=outmin'60 
c  Data  set  definitions 

read(61  ,*)(option(l),l=1  ,nops) 
c  Graphics,  column  data  to  sent  to  graphics  files. 

read(61,*)nslices 
read(61  ,*)(slice(K),K=1  .nslices) 

do  598 1=1, NC 
READ(61,*)NOCOP(l) 

NOCOPV=NOCOP(l) 

READ(61  ,*)(CELLOP(l,M),M=1  ,NOCOPV) 

598  CONTINUE 

do  601  1=1  ,NC 
M=1 

CELL  =  CELLOP(I.M) 

NNLMA=NNL(I) 

DO  602  J=1,NNLMA 
IF(CELLEQ.J)THEN 


VLl 


I 

DUMMY(I,J)=1 

M=M+1 

CELL=CELLOP(l,M) 

ELSE 

DUMMY(l,J)=0 

ENDIF 

602  CONTINUE 
601  CONTINUE 

c  read  time  in  minutes  at  which  output  is  to  start 
READ(61,')  OUTSTART 
OUTSTART=OUTSTART*60 


c 

c 

c 

c 


c 


350 

c 

c 

c 


c 


351 

c 

c 

c 

c 

c 

c 

c 

c 

c 


7 


8 


**«**«**«•«***«**«***«*******•**«#**•*•*«••»■*«•*•**•«*«*•**«***••* 


Mb.data 


Number  of  columns  for  which  mass  balance  data  is  to  be  written 
read{20,*)ncol 


Column  numbers  for  which  mass  balance  data  is  to  be  written 
do  350  l=1,ncol 
read(20,*)col{l) 
continue 

**t ***•**«*•«**« t****#***#**#*## #*««•***#****••*««••*•***•*•******* 

iso.vJata 

***•*«*•****«•*•»( . «*««**•.<****«*«***»* t «****«*•*#«*****••*««***** 

Number  of  columns  for  which  isotherm  is  to  be  identified 
read(21,*)isoncol 


Column  numbers  for  which  isotherm  is  to  be  identified 
do  351  l=1,isoncol 
read(21,*),isocol(l) 
continue 

*♦**♦***♦♦♦♦♦*♦***♦****♦***«**♦•**•*******•»******♦*•***********♦** 
INITIAL  SECTION 

♦****♦■*♦*♦*♦***#***♦***•*•****•**♦**♦♦**♦*♦**♦♦***********♦*****•** 

MILLINGTON  and  QUIRK  method  is  used  to  calculate 
hydraulic  conductivity  for  each  layer  from  a  given 
soil  moisture  characteristic  curve. 

NQJ=NQ1 

Layer  1 
D0  6I=1,NQJ 
IIJ=NQJ-I+1 
XII=X1(IIJ) 

TOPS=0.0 

BOTS=0.0 

D07J=1,NQJ 

JF=NQJ-J+1 

YJJ=Y1(JF) 

BOTS=((2*J-1)*YJJ**(-2))+BOTS 

11=1 

DO  8  J=II,NQJ 

JF=NQJ-J+1 

YJJ=Y1(JF) 

TOPS=((2*J+1-2*l)*YJJ**(-2))+TOPS 


Z1(JT)=STC0Nr(XII/ASR1)'T0PS/B0TS 

open(unit=2,file='zvar,status=‘unknown') 

6  write(2,‘)JT,Z1(JT),XII 

c  Layer  2 

NQJ=NQ2 
D0  9I=1.NQJ 
IIJ=NQJ-I+1 
XII=X2(IIJ) 

TOPS=0,0 

BOTS=0.0 

DO10J=1,NQJ 

JF=NQJ-J+1 

YJJ=Y2(JF) 

10  BOTS=((2‘J-1)*YJJ**(-2))+BOTS 
11=1 

DO  11  J=II,NQJ 

JF=NQJ-J+1 

YJJ=Y2(JF) 

1 1  TOPS=((2*J+1  -2*l)'YJJ**(-21)+TOPS 
JT=NQJ-I+1 

9  Z2(JT)=STCON2*(XII/ASR2)*TOPS/BOTS 

c  Calculate  gradients  on  suction/moisture  curves,  K/moisfure 

c  curve,  and  T/moisture  characteristic  curve. 

c  Calculate  average  incremental  gradients  of  suction-moisture 

c  and  K-moisture  curves, 

c  Layer  1 

D012I=1,.NQ1-1 

G1(I)=(Y1(I+1)-Y1(I))/(X1(I+1)-X1(I)) 

GZ1{I)=(Z1(I+1)-Z1(I))/(X1(I+1)-X1(1)) 

12  CONTINUE 

c  Layer  2 

D013I=1,NQ2-1 

G2(I)=(Y2(I+1  )-Y2(l))/(X2(l+1  )-X2(l)) 

GZ2{I)=(Z2(I+1  )-Z2(l))/(X2(l+1  )-X2(l)) 

13  CONTINUE 


c  T/0  curve. 

D014I=1,TF-1 

CVGR(I)=(WATECV(I+1)-WATECV(I))/(TEMCV(I+1)-TEMCV(I)) 
14  CONTINUE 

^  **************************«««**************»****************«*****«t* 

c  SET  PARAMETERS 

Q  ********************************************************************* 

c  Latent  heat  of  fusion,  units  J  g-1 

LA=333.2852 


c  DYNAMIC  SECTION 

Q  ******•♦******■***************♦**♦♦**♦•*♦*♦***•***♦*♦**♦*****♦*****♦** 

c  Set  number  of  iterations 

c  Iterations  =  no.  of  seconds  in  run/  iteration  period 


NOIT=NOSECR/IP 


read(49,*)NOTC,STAMBT,BTEMP 


DO  600  KKK=1,NOIT 
in=KKK 

DO  650 1=1, NC 
M=NNL(I) 

DO  651  J=1,M 
STATH(I,J)=ATHESP(I,J) 

STVATH(I,J)=VATHET(I,J) 

651  CONTINUE 

650  CONTINUE 

c  Second  marking  the  end  of  iteration 
NSEIT(KKK)=ITT*1P 

JJJ=0 

c  Set  up  loop  for  the  number  of  temp,  changes  in  data  set. 
DO  500  JJJ=1.NOTC 
if(JJJ.eq.1)then 
rewind  52 
rewind  48 
endif 

c  read  in  time  in  minutes  of  next  temperature  change 
read(52.*)TTMIN(JJJ) 
c  convert  this  value  to  seconds 
TTSEC(JJJ)=nMIN(JJJ)‘60 
read(48,*)AMBT(JJJ) 


if(TTSEC(JJJ).gt.NSEIT(KKK))then 

if(KKK.eq.1)then 

SEC(JJJ)=nSEC(JJJ) 

ATEMP(KKK)=STAMBT+(((AMBT(JJJ)-STAMBT)/SEC(JJJ)) 

&  *IP) 
endif 

if(KKK.gt.1)then 

ATEMP(KKK)=ATEMP(KKK-1  )+(((AMBT(JJJ)-ATEMP(KKK-1 )) 
&  /(TTSEC(JJJ)-NSEIT(KKK-1)))*IP) 
endif 

do  900  NN=1,nops 
if(option(NN).eq.31  )then 

open(unit=8,  file='atemp.data',  status='unknown') 
rewind  8 

open(unit=9,  file='Gatemp.data‘,  status='unknown') 
rewind  9 

write{8,*)'ATEMP=',ATEMP(KKK),KKK=',KKK 

write(9,*)ATEMP(KKK),KKK 

endif 

900  continue 


goto  700 
endif 


c  . 

if(nSEC(JJJ).eq.NSEIT(KKK))then 

ATEMP(KKK)=AMBT(JJJ) 

do  1000  NN=1,nops 

it(option(NN).eq.31  )then 

open(unit=8,file='atemp.data',status='unknown') 

rewind  3 

open(unit=9,  file='Gatemp.data',  status='unknown') 
rewind  9 

write(8,*)'ATEMP=',ATEMP(KKK);KKK=',KKK 

write(9.’)ATEMP(KKK),KKK 

endif 

1000  continue 
goto  700 
endif 

if(nSEC(JJJ).lt.NSEIT(KKK))goto500 


500  continue 

700  call  ku(NC,  LC1,  FC2.  KQ,  ATHESP, 

&  X1,X2,Z1,22,GZ1,GZ2,KP) 

call  kta(VATHET,  ATHESP,  CELVOL, 

&  sitp,  NNLMAX,  NC,  NNL,  ASR1 ,  ASR2,  LA, 

&  APKT,  KT,KI,NKI,KO,NKO,TKI,TKO,KS,OKG,IKG,SOILP) 

call  healrans(ATEMP,  CDEPTH,  slip,  sltpt2, 

&  VATHET,  NC,  NNL,  VHC,  LENGTH,  option, 

&  ATHESP,  NNLMAX,  BDT,  TEMCV,  WATECV, 

&  TF,  LA,  ORES1 ,  ORES2,  CVGR,  ATHET,  HFLUX,  NHFLUX, 

&  APKT,  STAMBT,  IP,  KKK,  LC1,  FC2,  ASRI,  ASR2, 

&  SHVATHET,CELVOL,ORIVOL,VADFAC,SOILP,SOIL1,TDEP,BTEMP,EXWAT) 

callmasstrans(NC,LC1,NQ1,ATHESP,X1,Y1,Gl,FC2,NNL, 

&  Y2G2,HT,KP,KQ,Z1,GZ1,Z2,GZ2,CDEPTH, LENGTH, 

&  sct,HPOT,AVK,FLUX,NFLUX,HYDG,IP,DEPTHP,TTCOM,dcp, 

&  BDT,  X2,  NQ2,  STCON1,  STCON2,  ASRI,  ASR2, 

&  EXSAT,  BWFLUX,  HTHETAC,  SHATHET,VATHET,COLHT, 

&  CDEP,EXP,vadfac,MAXCOLHT) 


c  call  massbal(ATHESP,  VATHET,  NC,  NNL,  STATH,  STVATH, 
c  &  NFLUX,  KKK,  EXSAT,  NSEIT,  BWFLUX,  HTHETAC, 
c  &  SHATHET,  SHVATHCT,  col,  ncol,  nops,  option) 


call  soli(NC,  NNL,  NNLMAX,  TTCOM,  set,  FS,  FDMAX,  FDMIN, 

&  FMAX,  FMIN,  FCUM,GAMMAS,GAMMAU,BETA,C,PHI,maxcolht, 

&  DEPTHP,HDISP,GELRAT,GDISP) 

if(mod(NSEIT(KKK),oufsec).eq.0.and.NSEIT(KKK).ge.OUTSTART) 

&  call  doutput(outsec,nops, 

&  option,NC,sltp,VATHET,ATHESP,sct,KP,NHFLUX,HFLUX, 


i 


&  VHC,FCUM,FMIN,FDMIN,FMAX,FDMAX.FS,APKT,KKK,NNLIVIAX, 

&  NOTC,  ATH1,  VATH1,  nSEC,  ATEMP,  BTFLUX,  KT,  nslices, 

&  slice, OUTTIM.DEPTHP.COLHT, EXP, SOILP,HDISP,EXWAT,GDISP) 


do  888  NN=1,nops 
if(option(NN).eq.35)then 

if(mod(NSEIT(KKK),outsec).eq.O)call  isotherm(NC,  sitp, 

&  TTCOM,dcp,KKK,NNL,ATEMP,ISO,outsec,isocol,isoncol, 
&  nops,option) 
endif 

888  continue 


600  continue 
DO  1=1, NC 

write(40,’)'total  downslope  displacement  at  surface' 
write(40,’)'column  ',l,  (GDISP(l)+HDISP(l))*100,'cm' 
enddo 
c 
c 
c 

stop 

end 


c  Subroutine  Cv.f 

c 

c  Calculates  volumetric  heat  capacity  of  a  soil  cell  given 
c  intrinsic  thermal  properties  of  soil  mineralogy,  and  relative 
c  proportions  of  soil  ice  and  water, 
c 

c  DEFINITIONS 

c  Volumetric  heat  capacity  (Cv)  is  the  amount  of  heat  required  to  change 
c  the  temperature  of  a  unit  volume  of  a  substance  by  1 K. 

c  Units  =  KJ  /  m-3 .  K 

c 

c  For  unfrozen  soil 

c  Cvu  =  d(Cs  +  CwW/100) 

c 

c  For  frozen  soil 

c  Cvf  =  df  (Cs  +  Cw  W/1 00  +  Ci  I/I  00) 

c 

c  Cvu  =  volumetric  heat  capacity  of  unfrozen  soil  (J  m-3  K-1 ) 
c  Cvf  =  volumetric  heat  capacity  of  frozen  soil  (m-3  K-1 ) 
c  d  =  dry  bulk  density  of  unfrozen  soil  (g  m-3) 

c  df  =  dry  bulk  density  of  frozen  soil  (g  m  -3) 

c  Cs  =  specific  heat  capacity  of  dry  mineral  matter  (J  g-1  K-1 ) 
c  Cw=  "  “  ■  ■  water  (J  g-1  K-1) 

c  Ci  =  “  ■  "  "  ice  (J  g  -1  K  -1) 

c  W  =  water  (volumetric  %)and  W/100  is  unit  %  ie.  ATHESP 

c  I  =  ice  (volumetric  %)  and  1/100  is  unit  %  ie.  VATHET 
c 

^  1 1  **«*«#*«**«****#*•*«*•*****•*  *^*«***«***«************************«»* 

subroutine  Cv(NC,NNL,VHC, ATHESP, VATHET.TT.NNLMAX) 

^  t***«*****«**«********«««*«**«**«»««*******«*************«********t«* 

c  DECLARATION  OF  VARIABLES 

^  *******«**•*•«*««**«******«********«****•*****««*«*****•**«********** 

DOUBLE  PRECISION  UDENS,  FDENS,  SHCSOIL,  VHC(5,60) 
DOUBLE  PRECISION  CW,  Cl,  ATHESP(5,60),VATHET(5,60) 

INTEGER  NC,  NNL(5),  H,  NNLMAX 

^  •********************t«***«****««***********************t***«******%* 

c  OPEN  STATEMENTS 

Q  #•**«****«***•***«**«****«**««*********«****«**********«***«**«****»* 

open(unit=59,  file='DATAIN/Cvin.data',  status='unknown') 
rewind  59 

^  t *««*****•**•«««***** t**«***«««*************«****«******************* 

c  READ  DATA  INPUT  FILE 

Q  ****************************t*****************************t********** 

READ(59,*)UDENS,  FDENS,  SHCSOIL,  CW,  Cl 
c  DYNAMIC  SECTION 

c  Values  of  ice  and  water  transferred  from  main  program  are  already  in 

c  %  volumetric  terms, 

c  Determine  whether  soil  cell  contains  ice. 

D01  l=.,NC 
M=NNL(1) 


Frozen  soil.  VHC  units  =  J  m  -3  K-1. 

IF(VATHET(l,J).GT.O.O)VHC(l,J)=FDENS* 

&  (SHCSOIL+(CW*ATHESP(l,J))+(CI*VATHET(l,J))) 


Unfrozen  soil.  VHC  units  =  J  m  -3  K-1. 

IF(VATHET(I,J).LE.O.O)VHC(I,J)=UDENS*(SHCSOIL+ 
&  (CW*ATHESP(I,J))) 


CONTINUE 

CONTINUE 


return 

end 


c 


Subroutine  freeze.f 


.  ****«•**•**«***«***«*«***«*****««*•«*********•***•***«*•*«•******«** 

c  Program  freeze.f 

................................................................ 

c  freeze  is  entered  if  the  net  heat  flux  for  a  cel!  is  negative 
c  as  calculated  in  heatrans.f 

subroutine  freeze(ATHESP,  TEMCV,  WATECV,  NC, 

&  TF,  sitp,  sltpt2,  CVGR,  VHC,  lA,  H.  LENGTH,  CDEPTH, 

&  VATHET,  BDT.  ORESI ,  ORES2,  ATHET,  option.  I,  J,  KKK, 

&  ASR1,ASR2,CELVOL,VADFAC,SOILP) 

c  DECLARATION  OF  VARIABLES 

^  •** t**********«**o*******««****««***************«***************«*** 

DOUBLE  PRECISION  ATHESP(5,60),  ATHET(5,60) 

DOUBLE  PRECISION  ATHESP1(5,60).  ATHET1(5,60) 

DOUBLE  PRECISION  VATHET(5.60),VATADD(5,60) 

DOUBLE  PRECISION  TEMCV(50),  WATECV{50) 

DOUBLE  PRECISION  sltpt2(5,60),  sltp(5,60),  CVGR(50) 

DOUBLE  PRECISION  deltaT(5,60).  deltaW(5.60),  cmdt(5,60) 
DOUBLE  PRECISION  LAO(5.60),  VHC(5,60),  U 
DOUBLE  PRECISION  CH(5,60).VADFAC(5,60),NEWVOL(5,60) 
DOUBLE  PRECISION  LENGTH,  CDEPTH(5,60),  BDT 
DOUBLE  PRECISION  ORESI  (50),  ORES2(50),  MORES(50) 
DOUBLE  PRECISION  CHPERC(5,60),  TEMPCH(5,60),SOILP(5,60) 
DOUBLE  PRECISION  ASR1,  ASR2,  ADFAC(5,60),CELVOL(5,60) 


INTEGER  NC,  TF,  TFMAX,  option,  I,  J,  KKK 
ADFAC(I,J)=1,0 

^  **«******«'*********«•«’**•*«***•********•*•«****«*«******«******«•««* 

c  DYNAMIC  SECTION 

Q  t******************************************************************* 

c  Freezing  conditions  in  cell? 

IF((sltpt2(l,J)-273.1 5),  LT.  1  .Oe- 1 0)GOTO  1 2 
IF(slfpt2(l,J).LE.273.15)GOTO  12 
IF((sltpt2(l,  J)-273. 1 5).GT.  1  .Oe-1 0)THEN 
sltp(l,J)=sltpt2(l,J) 

GOTO  2 
ENDIF 


c  New  soil  temperature  at  end  of  iteration  has  been  calculated  by 
c  heatrans.f,  (sltpt2),  using  the  0/T  characteristic  curve, 
c  Establish  corresponding  soil  moisture  content  for  temperature 
c  sltpt2. 

c  The  value  returned  from  the  0/T  curve  is  adjusted  by  the  volume 
c  change  factor  for  the  cell.  This  means  that  the  absolute  volume  of 
c  water  at  a  given  temperature  is  the  same  whilst  the  volumetric 
c  percentage  will  be  decreased  if  expansion  has  taken  place. 

c  Calculate  dt  deltaT,  for  iteration  period.  (K) 

12  IF(sltp(l,J).GT.273.15)deltaT(l,J)=273.15-sltpt2(l,J) 

IF(sltp(l,J).LE.273.15)deltaT(l,J)=sltp(l,J}-sltpt2(l,J) 
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lF(deltaT(l,J).EQ.0.0D-10)THEN 
GOTO  2 
ENDIF 

IF(sltpt2(l,J).GE.TEMCV(1))ATHET(l,J)=WATECV(1) 

TFMAX=TF 

iF(sltpt2(l,J).LE.TEMCV(TFMAX))ATHET(I.J)=WATECV(TFMAX) 

D04K=1,TF-1 

IF(sltpt2(I.J).LEJEMCV(K).AND.sltpt2(l,J).GT.TEMCV(K+1)) 

&  ATHET(l,J)=WATECV(K)+CVGR(K)*(sltpt2(I.J)-TEMCV(K)) 


4  CONTINUE 

ATHET(I,J)=ATHET(I,J)/VADFAC(I,J) 

c  - 

c  Calculate  dW  deltaW. 

c  Calculate  mean  residual  moisture  content  and  adjust  for  any  volume 

c  change  that  has  occurs  based  on  the  absolute  residual  moisture 

c  content  being  unchanged. 

MORES(l)=((ORES1(l)+ORES2(l))/2)/VADFAC(l,J) 

c  Balance  checks,  ATHESP  must  not  fall  below  residual  moisture  content 
c  adjusted  for  cell  volume.  The  assumption  is  made  that  the  residual 
c  moisture  content  for  a  cell  is  an  absolute  volume  and  not  a  unit 
c  percentage  so  it  is  adjusted  for  any  change  in  volume. 

IF(DABS(ATHESP(I,J)-MORES(I)).LT.1.0D-10)ATHESP(I,J)=MORES(I) 
c  . no  further  freezing  unless  water  migrates  to  cell. 

IF(ATHESP(l,J).LT.MORES(l))ATHESP(l,J)=MORES(l) 
c  . no  further  freezing  unless  water  migrates  to  cell. 

c  Convert  ATHESP  from  %volume  to  volume  m3. 

ATHESP1(l,J)=CELVOL(l,J)*ATHESP(l,J) 

c  Balance  checks ,  ATHET  must  not  fall  below  residual  moisture  content. 
IF(DABS(ATHET(I,J)-MORES(I)).LT.1.0D-10)ATHET(I.J)=MORES(I) 


IF(ATHET(l,J).LT,MORES(l))ATHET(l,J)=MORES(l) 

c  Convert  ATHET  from  %volume  to  volume  m3. 
ATHET1(l,J)=CELVOL(l,J)*ATHET(i,J) 

c  dO,  deltaW  m3 

deltaW(l,J)=ATHESP1  (l,J)-ATHET1  (I.J) 

Q  ******************************************************************** 

c  cmdt,  thermal  component  of  energy  budget.  J  (m-3  K-1 ). 
cmdt(l,J)=((VHC(l,J)*CELVOL(I.J))*deltaT(I.J)) 

c  LAO,  hydrological  component  of  energy  budget.  J  (m-3). 
UO(l,J)=(LA*1000000.0)*deltaW(l,J) 

^  *****«*•«*«•***«•*••**•«*«««*•*«««•• t******************************* 

c  PHASE  CHANGE  LOGIC 


c 

c 

c 


dQsLadO  dt-cmdt 
dt 


c  If  change  in  moisture  content  is  negative,  (given  the  predicted 
c  temperature  change  and  isothermal  nature  of  pase  change  in  previous 
c  iteration),  phase  change  cannot  occur.  Energy  is  therefore  not 
c  released  as  a  hydrological  response  to  transient  thermal  conditions 
c  and  thus  has  no  effect  on  the  temperature  of  the  cell. 

lF(deltaW(l,J).LE.O.O)THEN 
c  . no  further  change  in  0  or  I 


sltp(l,J)=sltpt2(l,J) 
GOTO  2 
ENDIF 


^  ****************tk*********«««************«******««**************»** 

c  LAO  =  cmdt 

IF(DABS(LAO(l,J)-cmdt(l,J)).LT.1.0D-10)THEN 
c  . isothermal  phase  change  occurs. 

c  Temperature 

c  . reset  temperature  to  that  before  heat  transfer,  i.e.  temp 

c  processes. 

C  lC6 

c  VATHET(I,J)=VATHET(I,J)+(ATHESP(I,J)-ATHET(I,J)) 

VATADD(I,J)=(ATHESP(I,J)-ATHET(I,J))*1.09 
VATHET(I,J)=VATHET(I,J)+VATADD(I,J) 

c  Moisture 

ATHESP(I,J)=ATHET(I,J) 

c  . where  ATHET  is  defined  by  the  characteristic  curve  for 

c  sltpt2. 

c  determine  if  new  volumetric  total  exceeds  unity 


NEWVOL(l,J)=VATHET(l,J)+ATHESP(l,J)+SOILP(l,J) 

IF(NEWVOL(I.J).GT.1.0)THEN 

ADFAC(i,j)=NEWVOL(l,J) 

VATHET(I,J)=VATHET(I,J)/ADFAC(I,J) 

ATHESP(I,J)=ATHESP(I,J)/ADFAC(I,J) 

SOILP(l,J)=SOILP(l,J)/ADFAC(l,J) 

CDEPTH(I,J)=CDEPTH(I,J)*ADFAC(I,J) 

CELVOL(l,J)=CELVOL(l,J)*ADFAC{l,J) 

endif 

GOTO  2 
ENDIF 

c  COMMENT 

c  Moisture  content  is  set  to  ATHET  on  characteristic  cunre  as  defined 
c  by  sltpt2,  however  after  phase  change  release  of  energy,  the  temp, 
c  returns  to  slip.  Therefore  at  the  beginning  of  the  next  time  step, 

c  soil  water  and  temperature  are  not  in  accordance  with  the 

c  characteristic  curve,  ie.  the  soil  moisture  content  is  less  than 
c  normally  expected  given  the  cell  temperature.  0  in  next  iteration 

c  issettoATHET->ATHESP. 

c  - 

c  LAO<cmdt 

c  NEW  IF  STATEMENT 
IF(deltaW(l,J).GE.O.O)THEN 
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c  IF(LAO(I.J).LT.cmdt(l,J))THEN 

c  . no  isothermal  phase  change,  cell  cools.  Energy  released  by 

c  cooling  of  cell  (thermal  component),  is  partially  compensated  by 
c  phase  change  energy  release.  Temperature  is  not  maintained  at 
c  273. 1 5  K,  rather  the  cell  cools  by  an  amount  equivalent  to  LAO  given 

c  the  volumetric  heat  capacity  of  the  soil  mixture, 

c  Temperature 

c  Calculate  the  temperature  change  expected  by  an  input  of  1J. 

TEMPCH(I,J)=1.0/VHC(I,J) 

c  T emperature  change  given  the  cell  volume. 

CHPERC(l,J)=TEMPCH(l,J)/CELVOL(l,J) 
CH(l,J)=LAO{l,J)*CHPERC(l,J) 

Sltp(l,J)=sltpt2(l,J)+CH(I.J) 

c  Isothermal  ??? 

if(sltp(l,  J).gt.273. 1 5)slfp(l,  J)=273. 1 5 

VATADD{I,J)=(ATHESP(I,J)-ATHET{I,J))*1.09 
VATHET(I,J)=VATHET(I,J)+VATADD(I,J) 
c  VATHET(I,J)=VATHET(I,J)+(ATHESP(I,J)-ATHET(I,J)) 

ATHESP(I,J)=ATHET(I,J) 


c  determine  if  new  volumetric  total  exceeds  unity 

NEWVOL(l,J)=VATHET(l,J)+ATHESP(l,J)+SOILP(l,J) 

if(NEWVOL(l,J).GT.1.0)THEN 

ADFAC(i,i)=NEWVOL(l,J) 

VATHET(I,J)=VATHET(I,J)/ADFAC(I,J) 

ATHESP(I,J)=ATHESP(I,J)/ADFAC(I,J) 

SOILP(l,J)=SOILP(l,J)/ADFAC(l,J) 

CDEPTH(I,J)=CDEPTH(I,J)*ADFAC(I,J) 

CELVOL(l,J)=CELVOL(l,J)*ADFAC(l,J) 

endif 

ENDIF 

^  **«***«*««««**«*•«**«******«****  ft**<r***ir****t**fttt*««**ii««««-k*li««-fc*«* 

2  return 

end 


XVlll 


c  Subroutine  heatrans.f 

c  Subroutine  calculates  flux  of  heat  energy  in  one  dimension  on  basis 
c  of  thermal  conductivity  of  medium  and  temperature  gradient, 
c  Temperature  of  cells  are  updated  given  the  volumetric  heat  capacity 
c  of  the  medium,  modified  on  each  iteration  for  changing  ice  and  water 
c  content, 

c 
c 

c  CdT  =  -d(kdT) 

c  dt  dx  dx 

c 

c  where 

c  C  =  volumetric  heat  capacity  J  m  -3  K  -1 
c  X  =  depth  m 

c  t  =  time  sec 

c  k  =  thermal  conductivity  Jsec  m  -1  K-t 

c  T  =  temperature  K 

c 

c  and 

c  1  W  =  1  J  sec 

c 
c 
c 

c  Subroutine  contains  time  element  for  thermal  conductivity 
c  Kt.  W  m-1  K-1  where  1 W  =  1J  sec 
c  i.e.  multiply  thermal  conductivity  entered  in  input  file 

c  by  iteration  period  in  seconds. 

-  ****•**•***•***«***«•*******•* ********************************* 

C 

subroutine  heatrans(ATEMP,  CDEPTH,  sitp,  sltpt2, 

&  VATHET,  NC,  NNL,  VHC,  LENGTH,  option, 

&  ATHESP,  NNLMAX,  BDT,  TEMCV,  WATECV,  TF,  LA. 

&  ORES1 ,  ORES2,  CVGR,  ATHET,  HFLUX,  NHFLUX,  APKT, 

&  STAMBT,IP,  KKK,  LC1 ,  FC2,  ASR1 ,  ASR2, 

&  SHVATHET,CELVOL,ORIVOL,VADFAC,SOILP.SOIL1,TDEP,BTEMP,EXWAT) 

^  *  «••■*•««««**«««****#*«****#******«'*****•**************************** 

c  DECLARATION  OF  VARIABLES 

^  ****«***********«***«*«*****t********««***************«***********«*« 

DOUBLE  PRECISION  ATEMP(1468800),  sltp(5,60),  sltpt2(5,60) 

DOUBLE  PRECISION  CDEPTH(5,60),  VATHET(5,60),  AVKT(5,60) 

DOUBLE  PRECISION  NHFLUX(5,60),  BTFLUX{50),  HFLUX(5,60) 

DOUBLE  PRECISION  LENGTH,  TCHA(5,60).CELVOL(5,60) 

DOUBLE  PRECISION  ATHESP(5,60),  ATHET(5,60).  BDT 
DOUBLE  PRECISION  VHC(5,60),  WAnS(5.60) 

DOUBLE  PRECISION  TEMCV(50).  WATECV(50).  LA,  CVGR(50),TDEP 
DOUBLE  PRECISION  ORES1(50),  ORES2(50).EXWAT(5,60) 

DOUBLE  PRECISION  APKT(5,60),  STAMBT.BTEMP 
DOUBLE  PRECISION  ASRI,  ASR2,SOILP(5,60),SOIL1(5,60) 

DOUBLE  PRECISION  SHVATHET(5,60),VADFAC(5.60),ORIVOL(5,60) 


c 


INTEGER  NC,  NNL(5),  option,  NNLMAX 
INTEGER  LC1(50),  FC2(50) 

INTEGER  TF,  l,J,IP,  KKK 


DYNAMIC  SECTION 


c 

c  Call  subroutine  Cv.f  for  calculation  of  volumetric  heat  capacities 
c  in  each  cell  of  matrix. 

callCv(NC,NNL,VHC.ATHESP,VATHET,n,NNLMAX) 

DO  666 1=1, NC 
M=  NNL(I) 

D0  667J=1,M 

SHVATHET(I.J)=VATHET(I,J) 

c  calculate  volume  adjustment  factor  (VA;:)FAC)  before  going  through 
c  heatrans  and  on  to  thaw  or  freeze 

VADFAC(l,J)=CELVOL(l,J)/ORIVOL(l,J) 

667  CONTINUE 
666  CONTINUE 

c  SURFACE  BOUNDARY  CONDITIONS 

c  Mean  thermal  conductivity  across  surface  interface. 

DO  1 1=1, NC 
J=1 

AVKT(I,J)=APKT{1,J) 

c  Surface  boundary  heat  flux, 

c  Divide  by  0.5  x  cell  depth. 

HFLUX(l,J)=((AVKT(!,J)*IP)*((ATEMP(KKK)-sltp(l,J))/ 

&  (0.5*CDEPTH(I,J))))*(LENGTH*BDT) 

1  CONTINUE 

c  PROFILE  HEAT  FLUX 

c  Mean  thermal  conductivities  across  cell  boundaries. 

DO  2  1=1, NC 
M=NNL(I) 

DO  3  J=2,M 

AVKT(I,J)=((APKT(I,J)*CDEPTH(I.J))+(APKT(I,J-1)* 

&  CDEPTH(I,J-1)))/(CDEPTH(I,J)+CDEPTH(I.J-1)) 


c  Heat  flux  across  cell  boundaries.  Jsec  cm-3  -K. 

HFLUX(l,J)=((AVKT(l,J)*IP)*((sltp(l,J-1)-sltp(l,J))/ 

&  ((CDEPTH(I,J)+CDEPTH(I,J-1))/2)))*(LENGTH*BDT) 

3  CONTINUE 

2  CONTINUE 

c  Basal  heat  flux.  Jsec  m-3  -K 

DO  4 1=1, NC 
J=NNL(1) 

c  btemp  set  as  input  basal  temperature 
c  also  possible  to  set  btflux  to  zero  if  do  not  want  basal  heat  supply 
btflux(i)=((avkt(i,j)*ip)’((sltp(i,j)-btemp)/cdepth(i,j)))* 

&  (length'bdt) 

4  CONTINUE 


c 


Net  heat  flux  across  all  cell  boundaries.  J  cm-3  -K. 


DO  5 1=1, NC 
M=NNUI) 
D0  6J=1,M 


c  Net  flux  of  heat  across  all  traudaries 

IF(J.NE.NNL(I))NHFLUX(I,J)=HFLUX(I,J)-HFLUX(I,J+1) 

c  Basal  cell  net  heat  flux. 

IF(J.EQ.NNL(I))NHFLUX(I.J)=HFLUX(I.J)-BTFLUX(I) 

6  CONTINUE 
5  CONTINUE 

c  T emperature  change  of  all  cells  in  matrix 

DO  7 1=1, NC 
M=NNL(I) 

D0  8J=1,M 

c  Calculate  Joules  required  to  change  temperature  by  1 K 
WAnS(l,J)=1.0/VHC(l,J) 

T^mnArstiirP  phannp 

TCHA(I,J)=(WATTS(I,J)/(LENGTH‘CDEPTH(I,J)*BDT))*NHFLUX(I,J) 

c  New  soil  temperature. 

sltpt2(l,J)=sltp(l,J)+TCHA(l,J) 

8  CONTINUE 

7  CONTINUE 

^  ****•***«****«* t*****«***«***««*»*«***«************«****«****«t******* 

208  DO  888  1=1, NC 

M=NNL(I) 

DO  999  J=1,M 

IF(NHFLUX(l,J).GT.0.0)THEN 
call  thaw(sltpt2,  TF,  LENGTH, 

&  CDEPTH,  BDT,  ATHESP,  VATHET,  VHC,  sitp,  NC,  LA, 

&  NHFLUX,  I,  J,  KKK,CELVOL,VADFAC,SOILP,ORIVOL,SOIL1 , 

&  ASR1,ASR2,LC1,FC2,TDEP,EXWAT) 

ENDIF 

c 

IF(NHFLUX(l,J).LE,0.0)THEN 

call  freeze{ATHESP,  TEMCV,  WATECV,  NC,  TF,  slip, 

&  sltpt2,  CVGR,  VHC,  LA,  TT,  LENGTH.  CDEPTH,  VATHET, 

&  BDT,  ORES1 ,  ORES2.  ATHET,  option,  I,  J,  KKK, 

&  ASR1,ASR2,CELVOL,VADFAC,SOILP) 

ENDIF 


999  CONTINUE 
888  CONTINUE 

Q  ********************************************************************** 

return 

end 


subroutine  isotherm.f 


This  subroutine  defines  and  writes  out  the  depth  of  the  273. 15K 
isotherm  to  a  file  which  may  then  be  convert^  into  a  graphical 
plot. 


subroutine  isotherm(NC,  sitp,  TTCOM,  dcp, 

KKK,  NNL,  ATEMP,ISO,outsec,isocol,isoncol, 
nops, option) 

DECLARATION  OF  VARIABLES 

DOUBLE  PRECISION  sltp(5,60),  ATEMP(1 468800) 
DOUBLE  PRECISION  TTCOM(5,60),  dcp{5.60) 
DOUBLE  PRECISION  ISO(5,60).  DISO(5.60) 

DOUBLE  PRECISION  SISO(5,60).TCH(5,60),TG(5.60) 
DOUBLE  PRECISION  TISO(525600).  BISO(525600) 


INTEGER  KKK,  NC,  NNL(5) 
INTEGER  outsec,  MM 
INTEGER  isoncol,  isocol(50) 
INTEGER  nops,  option(50) 


SURFACE  . 


!*#-•«««##•#*«*•*#•«••*•«•*•*•••••****•*•*•**«****** 


iLLS 


Does  273. 15K  isotherm  lie  between  surface  and  first  compn.  point? 
DO  1  1=1, isoncol 


M=NNL(I) 

D0  2J=1,M 

ll=isocol(l) 

if(J.eq.1)then 

if(sltp(ll,J).le.273.15,and.ATEMP(KKK).gt.273.15)then 

Temperature  gradient 
TG(ll,J)=(ATEMP(KKK)-sltp(ll,J))/dcp(ll,j) 

What  is  the  temperature  difference  between  surface  and  273.1 5K? 
TCH(II,J)=ATEMP(KKK)-273.15 

At  what  depth  from  the  surface  does  the  isotherm  occur? 
ISO(ll,J)=TCH(ll,J)n’G(ll,J) 

if(ISO(ll,J).le.O.O)goto  140 

Write  isotherm  data  to  files  with  suffix  300-399 

do  506  NN=1,nops 
if(option(NN).eq.35)then 
MM=(ll+300) 


XX 


write(MM,*)ISO(ll,J),KKK 

endif 

506  continue 
140  endif 


c  - 

if(sltp(l  I  ,J).gt.273. 1 5.and.ATEMP(KKK).le.273. 1 5)then 
TG(II.J)=(sltp(ll,J)-ATEMP(KKK))/dcp(ll.j) 
TCH(II,J)=273.15-ATEMP(KKK) 
SISO(II.J)=TCH(ll,J)/TG(ll,J) 

IS0(II.J)=S1S0(II,J) 

if(ISO(ll,J).le.0.0)goto150 


c  Write  isotherm  data  to  files  with  suffix  300-400 
do  501  NN=1,nops 
if(option(NN).eq.35)then 
MM=(I4300) 
write(MM,*)ISO(ll,J),KKK 
endif 

501  continue 

150  endif 

endif 


^  *****•«««****«****•*«**«*«***#*«*«*#*«******•*••****•**••***«•*** 

c  PROFILE  CELLS 


if(J.gt.1)then 

c  Does  273.1 5K  isotherm  lie  between  computation  points  of  cells 
c  J-1  and  J. 

if(sltp(ll,J).le.273.15,and.sltp(ll,J-1).gt.273.15)then 

c  Temperature  gradient 

TG(ll,J)=(sltp(ll,J-1)-sltp(ll,J))/(dcp(ll,J)+dcp(II.J-1)) 


c  What  is  the  temperature  difference  between  cell  J-1  and  273. 1 5? 
TCH(ll,J)=sltp(tl,J-1)-273.15 

c  At  what  depth  from  compn.  point  of  cell  J-1  does  273.1 5  isotherm 
c  occur? 

DISO(ll,J)=TCH(ll,J)/TG(ll,J) 

c  Depth  of  isotherm  from  surface. 

ISO(ll,J)=nCOM(ll,J-1)+DISO(ll,J) 

if(ISO(ll,J)  :e.0.0)goto  160 


do  502  NN=1,nops 


xxiii 


if(option(NN).eq.35)then 

MM=(ll+300) 

write(MM;)ISO(ll,J),KKK 

endif 

continue 


endif 


if(sltp(II.J).gt.273.15.and.sltp(ll,J-1).le.273.15)then 

TG(ll,J)=(sltp<ll,J)-sltp(ll,J-1))/(dcp(ll,J)+dcp(ll,J-1)) 

TCH(il.J)=273.15-sltp(ll,J-1) 

DISO(IU)=TCH(II.J)/TG(ll,J) 

lSO(ll,J)=TTCOM(ll,J)+DISO(ll,J) 

if(ISO(ll,J).le.O.O)goto  170 


do  503  NN=1,nops 
it(option(NN).eq.35)then 

MM=(ll-rC00) 

write(MM,*)ISO(ll,J),KKK 

endif 

continue 


endif 


endif 


continue 

continue 


return 

end 


c 


Subroutine  kta.f 


c  Sub  routine  kta.f  calculates  temperature  dependent  thermal 

c  conductivity  (W  m-1  K-1)  for  each  cell  of  profile.  The  system  is 

c  assumed  to  comprise  three  phases,  soil  solids,  water  and  ice. 


subroutine  kta(V.ATHET,  ATHESP.CELVOL, 

&  sitp,  NNLMAX,  NC,  NNL,  ASR1 ,  ASR2,  LA, 

&  APKT,KT,Kl,NKI,KO,NKO,TKI.TKO,KS,OKG,lKG,SOILP) 


c  DECLARATION  OF  VARIABLES 


DOUBLE  PRECISION  KI(20),TKI(20i,KO(20),TKO{20) 

DOUBLE  PRECISION  VATHET(5,6d),  ATHESP{5,60),CELVOL(5,60) 
DOUBLE  PRECISION  KS(5,60),  sltp(5,6u„  SOILV(5,60) 

DOUBLE  PRECISION  KOK(5,60),  KIK(5,60),  AGR1,  ASR2 
DOUBLE  PRECISION  VATH1(5,6G),  ATH1(5,60),SOILP(5,60) 
DOUBLE  PRECISION  SFI(5,60),  SFO,  SFS(5,60) 

DOUBLE  PRECISION  KT(5,60),  APKT(5,60),LA,IKG(20),OKG(20) 


INTEGER  NKI,  NKO,  NC,  NNL(5),  NNLMAX,  NKOMAX,  NKIMAX 

Q  t ***«**•«**#**««*««**•• t«««*«««*«««***«t*««**«»*****««*ft«««t*t *•««*«* 

c  DYNAMIC  SECTION 

^  *#**«*■•*»««•**«•»**•*****#*«*«*••««••*«*•«•**•**«•««*«**•***•»«•*•*«* 

c  Derive  volume  of  soii  soilds  per  cell 

DO  10 1=1, NC 
M=NNL(I) 

DO  20  J=1,M 

SOILV(l,J)=CELVOL(l,J)'SOILP(l,J) 

c  Derive  volume  of  ice  per  cell 

VATH1(l,J)=CELVOL(l,J)*VATHET(l,J) 

c  Derive  volume  of  water  per  cell 

ATH1(l,J)=CELVOL(l,J)'ATHESP(l,J) 

20  continue 
10  continue 


c  Extrapolation  from  curve 

c  Water 

DO  40 1=1, NC 
M=NNL(I) 

DO  41  J=1,M 

if(sltp(I.J).GE.TKO(1))KOK(l,J)=KO(1) 

NKOMAX=NKO 

if(sltp(l,J).LE.TKO(NKOMAX))KOK(l,J)=KO(NKOMAX) 

do4K=1,NKO-1 

if(sltp(l,J).LE.TKO(K).AND.sltp(l,J).GT.TKO(K+1)) 

&  KOK(l,J)=KO(K)+OKG(K)*(sltp(l,J)-TKO(K)) 

4  continue 

41  continue 

40  continue 


XXV 


c  Ice 

DO  60 1=1, NC 
M=NNL(l) 

DO70J=1,M 

if(sltp(l,J).GE.TKI(1))KIK(l,J)=K)(1) 

NKIMAX=NKI 

if(sltp(l,J).LE.TKI(NKIMAX))KIK(l,J)=Kl(NKlMAX) 


do5K=1,NKI-1 

if(sltp(I.J)-LE.TKI(K).AND.sltp(l,J).GT.TKI(K+1)) 

&  KIK(i,J)=KI(K)+IKG(K)*(sltp(l.vJ)-TKI(K)) 

5  continue 

70  continue 

60  continue 

c  Derive  shape  factor  for  water 

SFO=  1.0000 

c  Derive  shape  factor  for  ice 

DO  50 1=1, NC 
M=NNL(I) 

DO  51  J=1,M 
SF((I,J)=(1.0D00/3.0D00) 

&  *((1. 0000/(1. 0O00+((KIK(l,J)/KOK(l,J))-1.0D00)*0.125)) 

&  +(1  0D00/(1.0D00+((KIK(l,J)/KOK(l,J))-1.0D00)*0.125)) 

&  +(1  ODOO/i  1 .0D00+((KIK(l,J)/KOK(l,J))-1 .0D00)*0.750))) 

c  Derive  shape  factor  for  soil  solids 

SFS(I,J)=(1. 0000/3.0000) 

&  *{(1  0D00/(1.0D00+{(KS(i,J)/KOK(l,J))-1.0D00)*0.125)) 
&  +(1  0000/(1. 0D00+((KS(l,J)/KOK(l,J))-1.0D00)*0.125)) 
&  +(1 .0D00/(1  .ODOO+((KS(I,J)/KOK(I,J))-1 .0D00)*0.750))) 

c  Calculate  thermal  conductivity  of  the  cell  W  m-1  K-1 
KT(I,J)=((ATH1(I,J)*K0K(I,J)*SF0) 

&  +(VATH1(l,J)*KIK(l,J)*SFI(l,J)) 

&  +(SOILV(l,J)*KS(l,J)*SFS(l,J))) 

&  /((ATH1(l,J)*SFO)+(VATH1(l,J)‘SFI(l,J)) 

&  +(SOILV(l,J)*SFS(l,J))) 


APKT(I,J)=KT(I,J) 


51  continue 
50  continue 


return 

end 


XXVI 


c  Subroutine  ku.f 

c  Subroutine  ku.f  calculates  hydraulic  conductivity 

c  for  all  cells  in  the  profile 


subroutine  ku(NC,  LC1,  FC2,  KQ,  ATHESP, 

&  X1,X2,Z1,Z2,GZ1.GZ2,KP) 

c  Z1  Z2 . curve  gradients  on  suction/moisture  curve. 

DOUBLE  PRECISION  ATHESP(5,60),X1(50).X2(50) 

DOUBLE  PRECISION  Z1(50),Z2(50),GZ1(50),GZ2(50),KP(5,60) 

INTEGER  NC,  KQ,  LC1(50).  FC2(50) 

****«•*•***«******•*«*****«****••*****«*«*•***•*«*«*»«•*•**«***««*«* 
DYNAMIC  SECTION 

***♦*****•**♦************»*»****•»**••**»**»******••*•***********♦** 


c 

c 

c 


D01  1=1, NC 
M=LC1(I) 

D0  2J=1,M 

if(KQ.gt.20)write(428,*)KQ,athesp(i,j),i,j,kp(i,j) 

D0  3K=1,KQ-1 

IF(ATHESP(I,J).GT.X1(K).AND.ATHESP(I,J).LE.X1(K+1)) 
&  KP(I,J)=Z1(K)+GZ1{K)*(ATHESP(I,J)-X1(K)) 
IF(ATHESP(I,J).LE.X1  (1  ))KP(I,J)=Z1  (1 ) 
IF(ATHESP(I,J).GT.X1(20))KP(I,J)=Z1(20) 

3  CONTINUE 

2  CONTINUE 

1  CONTINUE 

DO  4 1=1, NC 
M=LC1(I) 

N=FC2{I) 

D0  5J=M.N 
D0  6K=1,KQ-1 

IF(ATHESP(I,J).GT.X2(K).AND.ATHESP(I,J).LE.X2(K+1)) 
&  KP(I,J)=Z2(K)+GZ2(K)*(ATHESP(I,J)-X2(K)) 
IF(ATHESP(I,J).LE.X2(1))KP(I,J)=Z2{1) 
IF(ATHESP(I,J).GT.X2(20))KP(I,J)=Z2(20) 

6  CONTINUE 

5  CONTINUE 

4  CONTINUE 


{20  Qo  Qo  Qo  Qo 


XXV 11 


c  Subroutine  masstrans.f 

c  Subroutine  masstrans.f,  is  used  to  reasses  the  distribution  of 
c  water  in  the  soil  profile  after  freezing  or  thawing  have  occurred, 
c  on  the  basis  of  Darcy's  Law. 

c 

subroutine  masstrans(NC,LC1,NQ1,ATHESP,Xl,Y1,G1,FC2,NNL, 
Y2,G2.HT,KP.KQ,Z1,GZ1,Z2,GZ2,CDEPTH,LENGTH, 
set, HPOT,AVK, FLUX, NFLUX,HYDG  'P  ""  TTCOM,dcp, 

BDT,  X2,  i\Q2,  STC0N1,  STC0N2,  a  EXSAT, 

BWFLUX,  HTHETAC,  SHATHET.VATHET,COLHT,CDEP,EXP, 
VADFAC,MAXCOLHT) 
c 
c 

DOUBLE  PRECISION  SHATHET(5,60).ASR1,ASR2,ATHESP(5,60) 
DOUBLE  PRECISION  ATHESPT(5.60).VATHET(5,60) 

DOUBLE  PRECISION  STC0N1,  STC0N2 
DOUBLE  PRECISION  Y1(50),Y2(50).G1(50),G2(50),X1(50),X2(50) 
DOUBLE  PRECISION  sct(5,60),  HPOT(5,60),CDEP(5,60),EXP(5,60) 
DOUBLE  PRECISION  HT(5,60),CDEPTH(5,60),LENGTH,BDT 
DOUBLE  PRECISION  AVK(5,60),KP(5,60),GZ1(50),GZ2(50) 
DOUBLE  PRECISION  FLUX(5,60),  NFLUX(5,60),BWFLUX(5) 
DOUBLE  PRECISION  COLHT(5),CUMHT,MAXCOLHT(5) 

DOUBLE  PRECISION  Z2(50),  Z1  (50).vadfac(5,60) 
c 

INTEGER  NC,  LC1(50),  FC2(50),  NQ1,  NQ2,  KQ,  NNL(5),IP 
c 

c  UNSATURATED  FLOW 

c  Soil  water  pressure  (suction) 

c  Layer  1 

DO  20I=1,NC 
M=LC1(I) 

DO  30J=1,M 

c  is  saturated  set  suction  to  zero 

IF(ATHESP(I,J).GE.ASR1)THEN 
SCT(l,J)=0.0 

ELSEIF(ATHESP(I,J).LT.X1(1))THEN 

SCT(I,J)=Y1(1) 

ELSE 
K=NQM 
DO  40  L=1,K 

IF(ATHESP(I,J).GE.X1  (L).AND.ATHESP(I.J).LT,X1  (L+ 1  ))THEN 
sct(l,J)=Yl  (L)+G1  (L)*(ATHESP(I,J)-X1  (L)) 

ENDIF 

40  CONTINUE 

ENDIF 

30  CONTINUE 

20  CONTINUE 


c  Layer  2 

DO  50 1=1, NC 
M=FC2(I) 

N=NNL(I) 

DO60J=M,N 

IF  (ATHESP(I,J).GE.ASR2)THEN 
SCT(l,J)=0.0 


xxviii 


ELSE1F(ATHESP(I.J).LT.X2(1))THEN 

SCT(I,J)=Y2(1) 

ELSE 

K=NQ2-1 

DO70L=1,K 

IF(ATHESP(I.J).GE.X2(L).AND.ATHESP(I,J).LT.X2(L+1))THEN 

sct(l,J)=Y2(L)+G2(L)*(ATHESP(l,J)-X2(L)) 

ENDIF 

70  CONTINUE 

ENDIF 

60  CONTINUE 

50  CONTINUE 


c  calculate  cell  mid-point  heights,  total  column  height 
c  and  volume  expansion  of  each  cell 
DO  151  1=1, NC 
M=NNL(I) 

CUMHT=0.0 
COLHT(I)=0.0 
D0 152  J=M,1.-1 

c  find  height  of  midpoint  of  cell 

CUMHT=CUMHT+(CDEPTH(l,J)/2.0) 
c  cell  expansion  calculated  as  value  in  mm 

EXP(I,J)=(CDEPTH(I,J)-CDEP(I.J))*1000 
HT(I.J)=CUMHT 

c  convert  CUMHT  to  height  of  column  so  far 

CUMHT=CUMHT+(CDEPTH(l,J)/2.0) 
COLHT(l)=COLHT(l)+CDEPTH(l,J) 

152  CONTINUE 

IF(COLHT(l).GT.MAXCOLHT(l))MAXCOLHT(l)=COLHT(l) 
151  CONTINUE 

c  Hydraulic  potential 

DO  80 1=1  ,NC 
DO  90J=1,NNL(I) 

HPOT(l,J)=sct(I.J)+HT(l,J) 

90  CONTINUE 

80  CONTINUE 

c 
c 

c  Hydraulic  conductivity 

c  Layer  1 

DO  100 1=1, NC 
M=LC1(I) 

DO110J=1,M 

IF(ATHESP(I,J).GE,ASR1)  THEN 
KP(l,J)=STCON1 

ELSEIF(ATHESP(I,J).LT,X1(1))THEN 

KP(I,J)=Z1(1) 

ELSE 

DO120K=1,KQ-1 

IF(ATHESP(I,J).GE.X1(K).AND.ATHESP(I,J).LT.X1{K+1)) 
&  KP(I  ,J)=Z1  (K)+GZ1  (K)*(ATHESP(I.J)-X1  (K)) 

120  CONTINUE 
ENDIF 

110  CONTINUE 
100  CONTINUE 


XXIX 


c  Layer  2 

DO  130 1=1, NC 
M=FC2(I) 

N=NNL(I) 

D0140J=M,N 

1F(ATHESP(I,J).GE.ASR2)THEN 

KP(l,J)=STCON2 

ELSE1F(ATHESP(I,J).LT.X2(1))THEN 

KP(I,J)=Z2(1) 

ELSE 

DO  150K=1,KQ-1 

IF(ATHESP(I,J).GE.X2(K).AND.ATHESP(I.J).LT.X2(K+1)) 
&  KP(I,J)=Z2(K)+GZ2(K)‘(ATHESP(I,J)-X2(K)) 

150  CONTINUE 
ENDIF 

140  CONTINUE 
130  CONTINUE 


c  Mean  hydraulic  conductivity  for  surface  cells 

DO  1401  l=1,NC 
J=1 

AVK(I,J)=KP(I,J) 

c  write(6,*)  'kp',KP(l,J),'avk',  AVK(1,J) 

1401  CONTINUE 


c  Mean  hydraulic  conductivity  for  profile  cells 
DO  1601=1,NC 
DO170J=2,NNL(l) 

c  if  freezing  has  started  use  upper  cell  conductivity  as  average 
c  conductivity  as  a  guide  to  ability  to  draw  up  water 
IF(VATHET(I,J-1).GT,0,0)THEN 
AVK(1,J)=KP(I  J-1) 

ELSE 

AVK(I,J)=((KP(I,J-1)*CDEPTH(I,J-1))+(KP(I,J)‘CDEPTH(I,J)))/ 
&  (CDEPTH(I,J-1)+CDEPTH(I,J)) 

ENDIF 

170  CONTINUE 

160  CONTINUE 

c 
c 

c  Convert  hydraulic  conductivity,  m  sec-1  into  m  lP-1 
DO  600 1=1, NC 
M=NNL(I) 

DO  601  J=1,M 
AVK(I,J)=AVK(I,J)*IP 
601  CONTINUE 

600  CONTINUE 

c  Calculate  flux  across  surface  boundary 
DO302l=1,NC 
J=1 

FLUX(l,J)=0.0 
302  CONTINUE 


XXX 


c  Calculate  Ilux  between  cells,  (vol  m-3  hr) 

D0 180 1=1, NC 
M=NNL(I) 

DO190J=2,M 

FLUX(l,J)=(((HPOT(I.J-1)-HPOT(I.J))*AVK(I.J))/((CDEPTH{l,J-1) 
&  +CDEPTH(l,J))/2))*(LENGTH*BOT) 
c  approach  to  stemming  downward  flux  if  cell  below  is 
c  saturated  with  water  or  water-fice. 

if(flux(i,j).gt.0.0.and.vadfac(i,j+1).gt.1.0)flux(i,j)=0.0 
190  CONTINUE 

180  CONTINUE 

c 
c 

c  Basal  water  flux 

DO  303 1=1, NC 
J=NNL(I) 

BWFLUX(I)=(AVK(I,J)*(0.5*CDEPTH(I.J)))*(LENGTH*BDT) 

303  CONTINUE 

c  Calculate  net  flux  for  each  cell 

DO  200 1=1, NC 
M=NNL(I) 

DO210J=1,M 

c  Profile  net  water  flux. 

if(J.NE.NNL(l))NFLUX(l,J)=FLUX(l,J)-FLUX(l,J+1) 
c  Basal  net  water  flux. 

if(J.EQ.NNL(l))NFLUX(l,J)=FLUX(l,J)-BWFLUX(l) 

210  CONTINJE 

200  CONTINUE 


c 

c  Calculate  new  soil  moisture  conditions 
DO  220 1=1, NC 
DO230J=1,NNU'l) 

ATHESPT(I,J)=ATHESP(I,J)+(NFLUX(I.J)/(CDEPTH(I,J)'LENGTH 
&  *BDT)) 

ATHESP(I,J)=ATHESPT(I,J) 

230  CONTINUE 
220  CONTINUE 
c 
c 

DO  240  1=1  .NC 
M=NNL(I) 

00  241  J=1,M 
SHATHET(I,J)=ATHESP(I,J) 

241  CONTINUE 
240  CONTINUE 


return 

end 


XXXI 


c  Subroutine  outpuLf 

c  Subroutine  OUTPUT.F  writes  data  to  files  according  to 

c  selection  of  option(l). 


subroutine  doutput(outsec,  nops,  option,  NC, 

&  sltp,  VATHET,  ATHESP.  set  KP,  NHFLUX,  HFLUX,  VHC, 

&  FCUM,  FMIN,  FDMIN,  FMAX,  FOMAX,  FS,  APKT,  KKK.NNLMAX, 

&  NOTC,  ATHI,  VATH1.  TTSEC.  ATEMP,  BTFLUX,  KT.nslices. 

&  slice, OUniM.DEPTHP.COLHT.EXP.SOILP.HDlSP.EXWAT.GDISP) 


c  DECLARATION  OF  VARIABLES 


DOUBLE  PRECISION  sttp(5,60),  VATHET(5,60),  ATHESP(5.60) 

DOUBLE  PRECISION  sct(5,60),  KP(5,60),  NHFLUX(5,60) 

DOUBLE  PRECISION  HFLUX(5,60),  VHC(5,60),SOILP(5,60) 

DOUBLE  PRECISION  FCUM(50),  FMIN(50),  FDMIN(50) 

DOUBLE  PRECISION  FMAX(50),  FDMAX(50),  FS(5,60) 

DOUBLE  PRECISION  APKr(5,60) 

DOUBLE  PRECISION  ATH1(5.60).  VATH1(5,60),  KT(5,60) 

DOUBLE  PRECISION  BTFLUX(50),DEPTHP(5),COLHT{5),EXP(5,60) 
DOUBLE  PRECISION  ATEMP(1468800),HDISP{5),GDISP(5),EXWAT(5.60) 

INTEGER  NC,  outsec,  nops,  option(50) 

INTEGER  KKK,  NNLMAX 
INTEGER  rTSEC{525600) 

INTEGER  nslices,  slice{50).OUniM 
COMMON  DUMMY(5,60) 


IF(OUTTIM.EQ.1)THEN 
do  100  NN=1,nops 
if(option(NN).eq.1)then 

open(unit=81 .  file='sltp.data'.s(alus=  unknown') 

rewind  81 

endif 


if(option(NN).eq.2)then 

open(unit=82,  file='ice.data‘,status='unknown') 

rewind  82 

endif 

if(option(NN).eq.3)then 

open(unit=83,  file='water.data',  status=‘unknown’) 

rewind  83 

endif 

if(option(NN}.eq.4)then 

open(unit=84,  file='nhflux.data',  status='unknown') 

rewind  84 

endif 

if(option(NN).eq.5)then 

opOT(unit=85,  file='hflux.data',  status='unknown') 

rewind  85 

endif 


XXX  ii 


if(option(NN).eq.6)then 

open(unit=86,  file='heatcap.data',  status=' unknown') 

rewind  86 

endif 

if(option(NN).eq.7)then 

open(unit=87,  file='thermcon.data‘,  status='unknown') 

rewind  87 

endif 

if(option(NN).eq.8)then 

open(unit=88,  file='hydcon.data',  sfafus='unknown') 

rewind  88 

endif 

if(option(NN).eq.9)then 

open(unit=89,  file='suction.data',  status='unknown') 

rewind  89 

endif 

if(option(NN).eq.  1 0)then 

open(unit=90,  fiie='fcum.data',  status='unknown') 

rewind  90 

endif 

if(option(NN).eq.11)fhen 

open(unit=91,  file='fmax.data',  status='unknown') 

rewind  91 

endif 

if(option(NN).eq.  1 2)then 

open(unit=92,  file='fmin.data',  status='unknown') 

rewind  92 

endif 

if(option(NN).eq.13)then 

open(unit=93,  file='fs,data',  status='unknown’) 

rewind  93 

endif 

if(option(NN).eq.  1 4)then 

open(unit=94,  file='ath1.data',  status='unknown') 

rewind  94 

endif 

if(option(NN).eq.  1 5)then 

open(unit=95,  file='vath1.data',  status='unknown') 

rewind  95 

endif 

if  (option(NN).eq.  1 6)then 

open(unit=96,  file='atemp.data',  status=’unknown') 

rewind  96 

endif 

if(option(NN).eq.17)then 

open(unit=97,file='btflux.data',  status='unknown') 


rewind  97 
endif 

if(option(NN).eq.20)then 

open(unit=70,  file='Gsltp.data‘,  status='unknown') 

rewind  70 

endif 

if(option(NN).eq.21)then 

open(!jnit=71,  file='Gice.data',  status='unknown') 

rewind  71 

endif 

if(option(NN).eq.22)then 

open(unit=72,  file=‘Gwater.data',  status='unknown') 

rewind  72 

endif 

if(option(NN).eq.23)then 

open(unit=73,  file='Gsuction.data',  status='unknown') 

rewind  73 

endif 

if(option(NN).eq.24)then 

open(unit=74,  file='Ghycond.data',  stafus='unknown') 

rewind  74 

endif 

if(option(NN).eq.25)then 

open(unit=75,  file='gtemp.data',  status='unknown') 

rewind  75 

endif 

if(option(NN).eq.26)then 

open(unit=76,  file='Gkt.data',status='unknown') 

rewind  76 

endif 

if(oplion(NN).eq.27)then 

open(unit=77,  file='Gapkt.data',status='unknown') 

rewind  77 

endif 

if(option(NN).eq.28)then 
open(unit=78,file='Gfs.data',status='unknown') 
rewind  78 
endif 

if(option(NN).eq.29)then 

open(unit=79,file='Gnhflux.data’,status='unknown') 

rewind  79 

endif 

if(option(NN).eq.30)then 

open(unit=80,file='Ghflux.data',status='unknown‘) 

rewind  80 

endif 


XXX IV 


if(option(NN).eq.33)then 
open(unit=33,file='Gexp.data‘,status='unknown') 
rewind  33 
endif 

if  ( option(NN).  eq.  34)then 

open(unit=34,file='Gcelexp.data',status='unknown') 

rewind  34 

endif 

if(option(NN).eq.36)then 

open(unit=36,file='Gsoilp.data',status='unknown') 

rewind  36 

endif 

if(option(NN).eq.37)then 

open(unit=37,file='Gdisplace.dala',slatus='unknown') 

rewind  37 

endif 

if(option(NN).eq.38)then 

open(unit=38,file='Gexwat.data',status='unknown') 

rewind  38 

endif 

100  continue 

OUTTlM=OUTTIM+1 

endif 

do  103  NN=1,NOPS 

c  OUTPUT  OPTIONS 

c  •**  SLTP  *** 

if(option(NN).eq.1)then 

wrile(81,11)KKK 

11  formal(//'Soil  temperature,  slip.  K',/, 

&  ’Iteration  =  ',i8) 

call  sieve(sltp,81,NC,NNLMAX) 

endif 

c  . 

c  **•  ICE  *** 

if(option(NN).eq.2)then 

write(82,13)KKK 

1 3  format(//'Soil  ice  content  m3  m-3',/, 

&  'Iteration  =  ',i8) 

call  sieve(VATHET,82,nc,nnlmax) 

endif 

c  - 

c  "*  WATER  *** 

if(option(NN).eq.3)then 


write(83,14)KKK 


XXXV 


1 4  format(//’Soil  moisture  content  m3  m-3',/, 

&  'Iteration  =  ',i8) 

call  sieve(ATHESP,83,NC,nnlmax) 

endif 

c  . 

c  *•*  NET  HEAT  FLUX  *** 

it(option(NN).eq.4)then 

write(84,16)KKK 

16  format(//'Net  heat  flux  W  m-3  K-1',/, 

&  'Iteration  =  ',i8) 

call  sieve(NHFLUX,84,nc,nnlmax) 

endif 

c  . 

c  ***  HEAT  FLUX 

if(option(NN).eq.5)fhen 

write(85,18)KKK 

1 8  format  (//'Heat  flux  W  m-3  K- 1 ',/, 

&  'Iteration  =  ',i8) 

call  sieve(HFLUX,85,nc,nnlmax) 

endif 

c  . 

c  '"VOLUMETRIC  HEAT  CAPACITY*** 

if(option(NN).eq.6)then 

write(86,20)KKK 

20  format(//'Volumetric  heat  capacity  J  m-3  K-1',/, 

&  'Iteration  =  ',i8) 


call  sieve(VHC,86,nc,nnlmax) 
endif 

c  . — . 

c  ***  THERMAL  CONDUCTIVITY  *** 

if(option(NN).eq.7)then 

write(87,22)KKK 

22  format(//Thermal  conductivity  W  m-3  sec-1  K-1 

&  'Iteration  =  ',i8) 

call  sieve(APKT,87,nc,nnlmax) 

endif 

c  - 

c  ***  HYDRAULIC  CONDUCTIVITY  *" 

if(option(NN).eq.8)then 

wiite(88,24)KKK 

24  formatl/yHydraulic conductivity  msec-1‘,/, 

&  'Iteration  =  ',i8) 


I  xxxvi 


call  sieve(KP,88,nc,nnlmax) 
endif 

c  . 

c  ***  SUCTION  '** 

if(option{NN).eq.9)t‘‘3n 

write(89,26)KKK 

26  format(//'Soil  suction  m',/, 

&  'Iteration  =  ',i8) 

call  sieve2(sct,89,nc,nnlmax) 

endif 

c  . 

c  CUMULATIVE  “F"  RATIO  *" 

if(option(NN).eq.10)then 

write(90,28)KKK 

28  format(//'Cumuiative  F  ratio  per  column',/, 

&  'Iteration  =  ',i8) 

write(90,*)(FCUM(l),|r1,NC) 

endif 

c  . 

c  MAXIMUM  "F"  RATIO  *** 

if(option(NN),eq.11)then 

write(91,29)KKK 

29  fQrmat(//'Maximum  F  ratio  per  column',/, 

&  ‘Iteration  =  ',i8) 

write(91,*)(FMAX(l),l=1,NC) 

write(91,30) 

30  format(//'Dspth  of  maximum  F  ratio  from  surface  m'//) 
write(91,*)(FDMAX(l),l=1,NC) 

endif 

c  - 

c  *•’  MINIMUM  "F"  RATIO  *** 
if(option(NN).eq.12)then 

write(92,31)KKK 

31  format(//'Minimum  F  ratio  per  column',/, 

&  'Iteration  =  ',i8) 

write(92,*)(FMIN(l),l=1,NC) 

write(92,32) 

32  format(//'Depth  of  minimum  F  ratio  from  surface  m'//) 
write(92,‘)(FDMIN(l),l=1,NC) 

endif 

c  - 

c  •**  FACTOR  OF  SAFETY  *** 

if(option(NN).eq.13)then 

write(93,33)KKK 

33  format(//'Factor  of  Safety ',/, 

&  'Iteration  =  ',i8) 


call  sieve(FS,93,nc,nnlfnax) 
eridif 

c  . - — 

c  ***  ATH1  *** 

if(option(NN).eq.  1 4)then 

write(94,200)KKK 

200  format(//'Soil  water  content,  volume  m3'./. 

&  'lteration=',  18) 

call  sieve(ATH1,94,nc,nnlmax) 
endif 

c  . 

c  *’*  VATH1  *’* 

if(option(NN).eq.15)then 

write(95,202)KKK 

202  format(//'Soil  ice  content,  volume  m3',/, 

&  'lteration=',  18) 

call  sieve(VATH1,95,nc,nnlmax) 
endif 

c  . 

c  ’•*  ATEMP  *** 

if(option(NN).eq,16)then 

write(96,204)KKK 

204  format(//‘ Ambient  temperature  change  and  time',/, 

&  'lteration=',i8) 

write(96,*)ATEMP(l),TTSEC(l) 

endif 

c  . 

c  *’*  BTFLUX  *** 

if(option(NN),eq.17)then 

write(97,206)KKK 

206  format('Basal  heat  flux  W  m-3  s- 1 ',/, 

&  'Iteration  =',  i8) 

write  (97,*)(BTFLUX(I),I=1,NC) 
endif 

^  ******-****««**********************«***«****************4*********** 

c  GRAPHICS 

^  **************♦#***♦♦*****♦♦♦♦*****#***♦*************•*#**♦****♦♦** 

c  ***  Gsltp  *** 

c  *******  change  this  for  correct  output  ******* 

do  800  K=1,nslices 
if(option(NN).eq.20)then 

c  do35J=1,NNLMAX 

c  l=slice(K) 

c  call  sieve(70,*)  KKK,  sltp(l,J) 

call  sieve(sltp,70,nc,nnlmax) 
c  call  sieve(112.*)  atemp(kkk) 

35  continue 


endif 

c  . 

c  Gice  *” 

if(option(NN).eq.21)then 

call  sieve(VATHET,71  ,nc,nnlmax) 
endif 

c  . 

c  Gwafer  *** 

if(option(NN).eq.22)then 

call  sieve(ATHESP,72,nc,nnlmax) 

endif 

c  ”*  Gsuction  **' 

if  (option(NNj .  jq.23)then 

call  sieve2(sct.73,nc,nnlmax) 
endif 

c  . — - 

c  **’  Ghycond 

if(option(NN).eq,24)then 

call  sieve2(KP  74,nc,nnlmax) 
endif 

c  . 

c  gATEMP  *** 

if(option(NN).eq.25)then 

do  220  l=1,NOTC 

c  call  sieve(75;)TTSEC(l),ATEMP(n 

220  continue 

endif 

c  . 

c  *'*  gKT  *** 

if(option(NN).eq.26)fhen 

call  sieve2(KT,76,nc,nnlmax) 
endif 

c  . 

c  ***  gAPKT  *** 

if(option(NN).eq.27)then 

do  251  J=1,NNLMAX 
l=slice(K) 

c  call  sieve(77,*)APKT(l,J) 

251  continue 
endif 

c  . 

c  •**  FS  ’** 

if(option(NN).eq.28)then 

do252J=1.NNLMAX 

l=slice(K) 

c  call  sieve(78,*)FS(l,J) 

252  continue 


XXX  ix 


endif 

c  . 

c  ’**  NHFLUX  *•* 

if(op(ion(NN).6q.29)then 
c  do253J=1,NNLMAX 

c  l=slice(K) 

call  sieve(NHFLUX,79,nc,nnlmax) 
c  253  continue 

endif 

c  - 

c  HFLUX  *** 

if(option(NN).eq.30)then 
do  254J=1,NNLMAX 
l=slice(K) 

c  call  sieve(80,MHFLUX(U) 

254  continue 
endif 

c  . 

c  ***  COLUMN  HEIGHT  (CM)  *** 
if(option(NN).eq.33)fhen 
do  255 1=1,NC 

write(33,*)COLHT(l)*100.0,(COLHT(l)-OEPTHP(l))'100.0 

255  continue 
endif 

c  . 

c  ***  CELL  EXPANSION  (MM)  “• 
if  (option(NN),eq.34)then 
call  sieve(exp,34,nc,nnlmax) 
endif 

c  . 

c  ”*  SOIL  VOLUMETRIC  PERCENTAGE  *’* 
if(option(nn).eq,36)then 
call  sieve(soilp,36,nc,nnlinax) 
endif 

c  . 

c  ***  SOIL  DOWNSLOPE  DISPLACEMENT  POTENTIAL  (cm)*** 
if(option{nn).eq.37)then 
do  256 1=1. NC 

write(37,*)HDISP(l)*100.,GDiSP(l)*100..(GDISP(l)+HDISP(l)) 

&  *100. 

256  continue 
endif 

c  . 

c  *•*  EXPELLED  EXCESS  WATER  VOLUMES  **** 
if(option(nn).eq.38)then 
call  sieve(exwat,38,nc,nnimax) 
endif 


c  *•*  ISOTHERM  ••• 
800  continue 
103  continue 


return 

end 


c  subroutine  sieve 

c  subroutine  sieve  sifts  the  desired  arrays  and  outputs  the  cells 
c  specified  in  the  Out.  data  file  thereby  enabling  selective  output 
c  rather  than  the  entire  column  to  be  printed. 


subroutine  sieve{array,  outid.nc.nnimax) 
integer  nc,nnlmax,outid,  counter 
COMMON  DUMMY(5,60) 
double  precision  buffer(60),  array(5,60) 

do  1013  i=1,nc 
counter=0 
do  2013  j=1,nnlmax 
if(DUMMY(i,j).eq.1)then 
counter=counter+1 
buffer(counter)=array(i,j) 
endif 

2013  continue 

write(outid,'(30f9.4,x)')(buffer(col),  col=1  .counter) 
1013  continue 
return 
end 


c  subroutine  sieve2 


subroutine  sieve2(array,  outid.nc.nnimax) 
integer  nc.nnimax.outid,  counter 
COMMON  DUMMY(5.60) 
double  precision  buffer(60).  array(5.60) 

do  1013  i=1.nc 
counter=0 
do  2013  j=1.nnlmax 
if(DUMMY(i.j).eq.1)then 
counter=counter+1 
buffer(counter)=array(i.j) 
endif 

2013  continue 

c  write(outid.’(30f9.4.x)')(butfer(col).  col=1. counter) 

write(outid.')(buffer(col).  col=1. counter) 

1013  continue 
return 
end 


xli 


c  Subroutine  soli.f 

c  Subroutine  soli.f  calculates  solifluction  movement  potential  for  column  cells 

c 

c 

subroutine  soli(NC,  NNL,  NNLMAX,  TTCOM,sct,FS,i-DMAX, 

&  FDMIN,  FMAX,  FMIN,  FCUM.GAMMAS.GAMMAU, BETA, C, PHI, MAXCOLHT 
&  .DEPTHP.HDISP.GELRAT.GDISP) 

c  DECLARATION  OF  VARIABLES 

DOUBLE  PRECISION  C(5,60),  BETA(5) 

DOUBLE  PRECISION  GAMMAS.  GAMMAU,GELRAT,GDISP(5) 

DOUBLE  PRECISION  sct(5,60),  PHI(5,60),  FS(5,60) 

DOUBLE  PRECISION  TTCOM(5,60),  FMAX(50) 

DOUBLE  PRECISION  FMIN(50),  SMAX(50) 

DOUBLE  PRECISION  FDMAX(50),  FDMIN(50),DEPTHP(5) 

DOUBLE  PRECISION  FCUM(50),MAXCOLHT(5),HDISP(5) 

INTEGER  NC,  NNL(5),  NNLMAX 

^  ***«*«*««•*****«*•«*««•«**•*«**«********•**•**•******»** 

c  calculate  maximum  do\«nslope  displacement  due  to  heave 

c  resettlement  using  maximum  displacement  and  slope  angle 

DO  5  1=1, NC 

HD1SP(I)=(MAXC0LHT(I)-DEPTHP(I))  *  TAND(BETA(I)) 
GDISP(I)=HDISP(I)*GELRAT 
5  CONTINUE 

Q  •«•«**•***«•**«#•*«*••*«***<«*«*«***•****•#******»****»»********** 

c  Set  FMAX  and  FMIN  limits 

DO  1 1=1, NC 
FMAX(l)=-20000 
FMIN(I)=20000 
SMAX(l)=-20000 
1  CONTINUE 

c  Calculation  of  factor  of  safety. 

DO  3 1=1, NC 
M=NNL(I) 

D0  4J=1.M 
IF(sct(l,J).GE.O.O)THEN 

FS(l,J)=(C(l,J)+(GAMMAS*nCOM(I.J)*(DCOS(BETA(l)))"2 
&  -sct(l,J)*9.8)*DTAN(PHI(l,J)))/(GAMMAS*TTCOM(l,J) 

&  •DSIN(BETA(l))*DCOS(BETA(l))) 

ENDIF 

IF(sct(l,J).LT.O.O)THEN 

FS(!,J)=(C(l,J)+(GAMMAU*nCOM(I.J)*(DCOS(BETA(l)))**2 
&  -sct(l,J)*9.8)*DTAN(PHI(l,J)))/(GAMMAU*nCOM(l,J) 

&  •DSIN(BETA(l))*DCOS(BETA(l))) 

ENDIF 

IF(FS(I.J).GT.FMAX(I))FMAX(I)=FS(I.J) 

FDMAX(J)=TTCOM(I.J) 

IF(FS(I,J).LT.FMIN(I))FMIN(I)=FS(I,J) 

FDMIN(l)=nCOM(I.J) 


xlii 


lRp:Hl.J).GT.SMAX(l))SMAX(l)=sct(l,J) 

FCUM(I)=FCUM(I)+FS(I,J) 


4  CONTINUE 

3  CONTINUE 


return 

end 


xliii 


Subroutine  thaw.f 

Subroutine  thaw.f  calculates  the  partition  of  energy  between 
two  processes;- 

(i)  phase  change 

(ii)  thermal  change 

when  NHFLUX  is  +ve  and  thaw  of  soil  moisture 
is  a  likely  result. 


subroutine  thaw(sltpt2,  TF,  LENGTH, 

&  CDEPTH,  BDT,  ATHESP,  VATHET,  VHC,  sitp,  NC,  LA, 

&  NHFLUX,  I,  J,KKK,CELVOL,VADFAC,SOILP,ORIVOL,SOIL1  ,ASR1  ,ASR2, 
&  LC1,FC2,TDEP,EXWAT) 

t**#*********************************************************-****** 

DECURATION  OF  VARIABLES 

DOUBLE  PRECISION  sltpt2(5,60),  sltp(5,60),ASR1,ASR2 
DOUBLE  PRECISION  ATHESP(5,60),  ATHESP1(5,60) 

DOUBLE  PRECISION  VATHET(5,60),  VATH1(5,60) 

DOUBLE  PRECISION  LENGTH,  CDEPTH(5,60),  BDT 
DOUBLE  PRECISION  deltaW(5,60),  deitaT(5,60),Et(5,60) 

DOUBLE  PRECISION  LA,  VHC(5,60),CELVOL(5,60),SOIL1(5,60) 

DOUBLE  PRECISION  SOILP(5,60),volch(5,60),ORIVOL(5,60) 

DOUBLE  PRECISION  VADFAC(5,60),TDEP 
DOUBLE  PRECISION  EXWAT(5,60),  WAT(5,60) 

INTEGER  LC1(50),FC2(50),TF,  NC,  I,  J,  KKK 

•***'t*t****»****«******************************************t*t***** 


Is  energy  partition  between  phase  change  and  thermal  heating 
required,  i.e.  is  ice  present. 

IF(VATHET(l,J).LE.0.0D-06.OR.SLTPT2(l,J).LT.TDEP)THEN 

sltp(l,J)=sltpt2(I.J) 

ELSEIF(SLTPT2(I,J).GE.TDEP.AND.VATHET(I,J).GT.0.0)THEN 

sltp(i,j)=tdep 

««««*«*«ir«’*#*«**#******************»*****«************************ 

Change  in  temperature  for  iteration  period  predicted  by  heatrans.f 
this  is  heat  energy  that  will  be  used  to  melt  ice 
deltaT(l,J)=sltpt2(I.J)-TDEP 

calculate  thermal  energy  available  for  melting  ice 
Et(I.J)  =  ((VHC(l,J)*CELVOL(l,J))*DeltaT(l,J)) 


HYDROLOGY 

Convert  Athesp  and  Athet  from  %  volume  to  volume  m3 
ATHESP1(l,J)=CELVOL(l,J)*ATHESP(l,J) 

Volumetric  change  in  moisture  content  -  ice  melted  (wat  equiv)  m3 
deltaW(l,  J)=Et(i,j)/(U*  1 000000.0) 

Convert  VATHET  ice  content  %  info  a  volume  m3 
VATH1(l,J)=CELVOL(l,J)*VATHET(l,J) 


check  that  ice  is  available  for  melting 
IF(VATH1(l,J).GT,(deltaw(I.J)*1.09))THEN 
VATH1(l,J)=VATH1(l,J)-(deltaW(l,J)*1.09) 
ELSE 

DELTAW(I,J)=VATH1(I.J)/1.09 

VATH1(I.J)=0.0 

ENDIF 

ATHESP1  (l,J)=ATHESP1  (l,J)+DELTAW(l,J) 

ATHESP(I.J)=ATHESP1(l,J)/CELVOL(l,J) 

VATHET(I,J)=VATH1(I,J)/CELV0L(I,J) 


calculate  if  any  volume  adjustment  is  required 
IF(VADFAC(I,J).GT.1.0)THEN 
IF(J.LE.LC1(l))WAT(l,J)=ASRrORiVOL(l,J) 
IF(J.GE.FC2(l))WAT(l,J)=ASR2*ORIVOL(l,J) 
include  lines  to  expel  water  content  that  is  above  the  original 
porosity  of  the  soil 

IF(ATHESP1(l,J).GT.WAT(l,J))THtN 

EXWAT(I,J)=EXWAT(1.J)+((ATHESP1(I,J)-WAT(I,J)) 

&  /ORIVOL(l,J)) 

ATHESP1(I,J)=WAT(I,J) 

ENDIF 

VOLCH(l,J)=CELVOL(l,JHATHESP1(l,J)+VA''H1(l,J)+SOIL1(l,J)) 

CELVOL(l,J)=CELVOL(l,J)-VOLCH(I.J) 

VADFAC(U)=CELVOL(l,J)/ORIVOL(l,J) 

IF(VADFAC(I,J).LT.1.0)THEN 

VADFAC(I,J)=1,0 

VOLCH(l,J)=CELVOL(l,J)-ORIVOL(l,J) 

CELVOL(l  J)=ORIVOL(l,J) 

ENDIF 

ATHESP(I,J)=ATHESP1(I,J)/CELV0L(I,J) 

VATHET(I.J)=VATH1(l,J)/CELVOL(I.J) 

S0ILP(I,J):-‘50IL1(I,J)/CELV0L(I,J) 

CDEPTH(l,J)=CDEPTH(l,J)-(VOLCH(l.J)/(BDrLENGTH)) 

endif 

else 

write(6,*)'slipped  out  in  thaw' 

ENDIF 


return 

end 


